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1  SUMMARY 


This  report  summarizes  the  activities  conducted  at  the  Center  for  Advanced  Power  and  Energy 
Research  (CAPEC)  from  2010  to  2014.  The  Center  was  fonned  to  serve  the  needs  of  the  Air 
Force  Research  Laboratory  (AFRL)  in  research  of  development  of  power  and  energy.  CAPEC 
was  a  collaborative  activity  between  AFRL  and  Wright  State  University  structured  through  a 
cooperative  research  agreement.  The  organizational  focuses  include: 

1 .  Modeling  of  plasma  physics 

2.  Modeling  fuel  cells 

3.  Testing  new  innovation  and  ideas  for  advanced  fuel  cells 

4.  Development  of  energy  related  issue  for  micro  air  vehicles  (MAVs). 

In  Section  2,  the  final  summary  of  present  work  on  plasma  modeling  is  discussed.  The  fuel  cell 
modeling  approach  is  given  in  Section  3.  The  development  in  fabrication  and  testing  of  advanced 
solid  oxide  fuel  cells  is  provided  in  Section  4.  The  characterization  of  MAV  propulsion  system 
will  be  given  in  Section  5. 


1 

Approved  for  public  release;  distribution  unlimited 


2  Plasma  Modeling 

2.1  Introduction 

Sustained  research  interest  and  achievements  in  flow  control  by  aerodynamics-electromagnetics 
interactions  have  grown  in  the  past  few  years.  Numerous  innovative  techniques  have  been 
developed  in  a  wide  range  of  applications  from  the  remote  energy  deposition,  electrical-thennal 
energy  conversion,  to  surface  plasma  generation  for  flow  control.  However,  the  surface  plasma 
actuators  are  the  most  frequently  adopted  technique  for  its  simplicity,  nonintrusive 
implementation  and  control  effectiveness  [1-5].  The  interaction  of  aerodynamics- 
electromagnetics  is  derived  from  the  three  basic  electromagnetic  properties  [6,7]:  The 
electrostatic  force  by  the  free-space  charge  separation  in  the  plasma  sheath  which  is  the 
cornerstone  of  dielectric  barrier  discharge  (DBD)  operation.  A  series  of  applications  devised  by 
Corke  et  al.  [1-2],  as  well  as,  Moreau  and  his  colleagues  [3]  have  been  successfully  demonstrated 
for  flow  control  at  subsonic  and  transonic  flow  regimes.  The  Joule  heating  occurs  for  all 
electrical  discharge  but  it’s  the  dominant  effect  for  direct  current  discharge  (DCD).  A  glow 
discharge  at  a  low  ambient  density  becomes  Corona  discharge  at  the  elevated  ambient  pressure 
condition.  The  thermal  plasma  actuator  is  best  suited  for  high  altitude  flight  and  closely 
associated  with  hypersonic  flows  [4,5],  The  third  mechanism  for  plasma  actuator  is  the  Lorentz 
acceleration  which  is  a  cross  product  of  an  externally  applied  magnetic  field  and  the  discharge 
current. 

In  surface  plasma  generation  by  the  electron  collision  process,  the  Townsend’s  mechanism 
controls  the  secondary  emission,  multiple  primary  avalanches,  and  ultimately  maintains  the 
discharge  [7,8],  For  this  reason,  the  classic  Townsend’s  similarity  law  for  electron  impact 
ionization  is  still  a  viable  fonnulation  of  the  complex  chemical-physical  process  to  simulate  the 
electrical  field  dominated  phenomena.  The  charged  particle  number  density  is  generally  limited 
to  an  order  of  magnitude  of  10 12  per  cubic  centimeter.  The  generated  plasma  consists  of  electrons 
in  a  highly  excited  state  but  the  heavy  ions  retain  the  thennodynamic  condition  of  their 
surrounding  environment.  Therefore,  the  weakly  ionized  gas  is  usually  far  from  thennodynamic 
equilibrium.  Meanwhile  the  drift  motion  of  charged  particles  and  diffusion,  including  the 
ambipolar  diffusion,  profoundly  modifies  the  transport  properties  of  the  ionized  medium. 

In  contrast  to  DCD,  the  DBD  is  maintained  by  an  alternating  electric  cunent.  The  DBD  operates 
on  a  self-limiting  process  through  the  reduced  electric  field  potential  by  the  surface  charge 
accumulation,  thus  prevents  the  corona-to-spark  transition  [7,8].  Specifically,  the  earlier 
outstanding  effort  by  Elisson  and  Kogelschlatz  [9]  has  identified  that  the  discharge  consists  of 
two  distinct  positive  Corona  streamers  and  diffusion  modes.  Enloe  et  al.  [10,1 1]  have  reaffirmed 
that  when  the  exposed  electrode  is  positively  biased  in  the  AC  cycle,  the  discharge  is 
characterized  by  a  streamer  like  structure.  These  microdischarges  have  a  short  life  span  of  about 
a  few  nanoseconds  and  with  the  random  temporal  and  spatial  structures.  In  the  negatively  biased 
phase  of  the  exposed  electrode,  it  acts  as  the  cathode  and  the  discharge  appears  as  a  more 
diffusive  structure  [9-11].  The  discharge  pattern  over  the  dielectric  surface  depends  on  the 
polarity  and  intensity  of  the  applied  electric  field,  as  well  as,  the  electric  pennittivity  of  the 
dielectrics  [7].  In  essence,  the  propagation  of  charged  particles  is  in  a  locked-step  to  the 
frequency  of  the  AC  field.  Meanwhile  the  induced  electrostatic  force  by  the  free-space  charge 
separation  during  the  AC  cycles  becomes  a  periodic  dynamic  event.  Nevertheless,  the  discharge 
phases  still  can  be  identified  as  avalanche,  streamer  formation,  a  subsequent  glow  discharge  and 
finally  quenching  of  the  microdischarge  on  the  electrodes  [7,8].  A  complex  physical 
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phenomenon  of  DBD  emerges;  a  wide  range  of  discharge  patterns  are  observed  depending  on  the 
gas  mixture  composition,  pressure,  electrodes  arrangement,  and  other  parameters.  However,  the 
global  structure  always  consists  of  the  continuous  diffusive  and  random  distributed  pulsing 
microdischarges  in  a  short  duration. 

An  enormous  amount  of  energy  is  needed  to  generate  localized  volumetric  plasma  that  must  have 
a  sufficient  charged  particle  number  density  for  strong  magneto-aerodynamic  interactions  [1-5]. 
For  examples,  the  ionization  potential  is  34  eV  for  electron  beam  [12],  65.7  eV  for  DCD,  and  81 
eV  per  ion-electron  pair  for  discharge  at  the  radio  frequency  [7,8].  The  ionization  potential 
always  underestimates  the  energy  requirement  in  applications,  because  the  nonequilibrium 
energy  cascades  to  vibration  excitation,  recombination,  and  attachment  processes.  In  the  electron 
impact  processes  for  ionization,  the  positive  and  negative  charged  ions  still  retain  their  ambient 
condition.  For  this  reason,  the  partially  ionized  is  often  identified  as  the  low-temperature  plasma 
with  a  charge  number  density  generally  limited  to  the  order  of  magnitude  of  10  '/cm  .  As  a 
weakly  partial  ionized  plasma,  the  electromagnetic  force  usually  exerts  a  small  perturbation  to 
the  mainstream  flow  and  the  thermodynamic  behavior  is  significantly  different  from  the  plasma 
generated  by  thennal  excitation  [7,1 1].  Therefore  the  plasma  actuator  for  flow  control  is  the  most 
effective  at  the  flow  bifurcations  such  as  the  onset  of  dynamic  stall,  laminar-turbulent  transition, 
vortical  separation  [1,2].  However  the  electromagnetic  effect  can  also  be  amplified  by  an 
externally  applied  magnetic  field  or  by  inviscid-viscous  interaction  at  the  leading  edge  of 
hypersonic  control  surfaces  [4,5]. 

The  nonequilibrium  chemical  kinetics  associated  with  the  DBD  in  atmosphere  is  well  known 
because  it  had  been  applied  as  an  ozone  generator  for  years.  Elisson  et  al.  [9]  have  identified 
plasma  chemistry  in  microdischarge  by  30  species  through  143  elementary  reactions.  In  a  more 
recent  work  by  Bogdanov  et  al.  [13],  the  chemical-physics  kinetics  of  atmospheric  plasma  have 
been  investigated  by  576  chemical  reactions  involving  vibrational  excitations  of  nitrogen  and 
oxygen,  ozone,  positive  and  negative  ions,  besides  oxide-nitrides.  The  complexity  of  the  internal 
degrees  of  excitations  includes  molecular  nitrogen  and  oxygen  at  few  quanta  above  ground  states; 
the  atomic  nitrogen,  (45, 2D,  2P) ,  oxygen  (3  p.'  s'  D) ,  the  charged  nitrogen  molecules 
(A3-z+u,B3ng,A'ng,c3nu)  and  oxy  gen  (x3i.+u,  a' a,  b'ti),  ozone  molecules  o3,  as  well  as,  negatively 
charged  ions  (n2o  ,no  , O', and  positively  charged  ions  (n+ ,o+ ,no+ ,no+,..) 
respectively.  Bogdanov  et  al  also  have  correctly  pointed  out  that  the  concentrations  of  the 
negatively  charged  ions  is  relatively  lower  than  the  positively  charged  counterpart,  but  the  effect 
of  the  negatively  charged  ions  can  still  exert  a  critical  directional  demarcation  for  the 

electrostatic  force  during  an  AC  cycle.  The  finding  by  Kim  et  al  [14]  has  also  shown  the  02 
plays  a  dominant  role  in  DBD  plasma  actuator.  Again,  DBD  is  a  nonequilibrium  transient 
discharge.  The  mean  energy  state  of  the  electrons  and  heavy  particles  is  considerably  different. 
For  this  reason  the  energy  of  charged  particles  can  be  effectively  transferred  to  the  internal  states 
of  molecule  or  atom. 

In  DBD  applications,  the  interaction  of  charged  particles  and  the  dielectrics  is  important. 
According  to  Golubovskii  et  al.  [15],  the  photoemission  is  negligibly  small  for  nitrogen  because 
the  energy  of  photons  in  the  visible  range  is  too  small  to  cause  the  emission  of  electrons.  Raizer 
[7]  also  has  pointed  out  that  the  photoemission  cannot  compete  with  electron  impact  ionization 
because  the  collision  cross  sections  of  molecular  and  atomic  oxygen  and  nitrogen  close  to  the 
emission  threshold  are  rather  high  (10‘  cm").  In  gist,  the  gas  in  DBD  occupies  only  few  quanta 


3 

Approved  for  public  release;  distribution  unlimited 


above  the  ground  state  that  is  capable  of  photoemission.  The  process  however  can  supply  the 
seeding  electrons  that  need  to  start  the  avalanche  in  streamer  propagation  [7.16]. 

The  other  important  processes  are  desorption  of  electrons  and  surface  recombination.  The 
electrons  can  be  desorbed  mostly  by  energy  cascading  to  vibrational  excitation  or  by  depletion  of 

the  metastable  molecule  N2(d3E+) .  However,  the  secondary  electrons  emission  of  electrons 
from  the  surface  during  an  AC  cycle  is  the  dominant  mechanism  and  leads  to  the  charge 
separation  in  the  cathode  layer.  Meanwhile,  the  surface  charge  accumulation  also  alters  the 
electric  potential  within  the  discharge  domain.  In  view  of  the  relatively  short  duration  of 
microdischarge  versus  the  microsecond  time  scale  of  the  AC  cycle;  the  instantaneous  surface 
recombination  is  a  plausible  and  acceptable  approximation.  In  this  connection,  one  needs  to  be 
mindful  of  the  fact  that  the  positive  charge  builds  up  on  the  dielectrics  far  downstream  of  the 
DBD  has  also  been  observed  by  Opaits  et  al.  [17].  They  also  noted  that  the  surface  potential 
coincided  with  the  biased  polarity  and  depletes  slowly  in  duration  of  tens  of  minutes  to  degrade 
the  performance  of  DBD.  Finally  and  most  importantly,  the  ionization  of  DBD  is  still  a  basic 
Townsend  process  [7,8].  This  observation  becomes  the  foundation  to  all  known  fonnulation  for 
ionization  by  electron  collisions. 

The  periodic  electrostatic  force  in  the  charge  separation  region  is  proportional  to  the  product  of 
the  net  balance  of  electrically  charged  number  density  and  electric  field  intensity  within  the 
plasma.  The  fonner  is  the  difference  between  the  number  densities  of  positively  and  negatively 

charged  ions  with  electrons,  pe  =  e(n (+  -  ne  -  ni  ) .  The  compatible  electric  field  intensity  is  also 
modified  by  the  distributions  of  the  charged  particle  density  or  the  charged  species  concentration 
[7].  On  the  other  hand,  the  distinct  differences  in  time  scales  in  the  fast  ionization  (10‘ns),  the 
slower  diffusion  phenomena  (10‘4s)  and  the  AC  cycle  (10‘3s)  are  the  key  elements  for  DBD 
operation  [13-18],  It  is  recognized  that  the  range  of  time  scales  varies  greatly  through  the  values 
from  tens  of  nanoseconds  to  a  few  millisecond  [16].  The  required  temporal  resolution  for 
experimental  measurements  and  computational  simulations  impose  a  significant  challenge 
[3,1 1,15-16].  In  computational  simulation,  the  spatial  resolution  requirement  for  multiple 
microdischarges  is  also  beyond  the  current  computational  capability.  Therefore  an  alternative 
approach  to  evaluate  the  low-temperature  weakly  ionized  gas  becomes  necessary. 

The  objectives  of  the  present  investigation  are  hopefully  to  clarify  some  of  the  conflicting 
explanations  by  putting  the  models  of  DCD  and  BDB  for  plasma  actuation  on  a  theoretical  frame 
based  on  the  drift-diffusion  and  chemical-physics  approximations.  Although  there  are 
deficiencies  in  ionization  by  electron  impact  process,  transport  property  by  the  drift-diffusion 
formulation,  and  bypassing  the  detailed  multiple  microdischarges  simulation,  but  the 
fundamental  physics  is  rigorously  described  by  the  electromagnetic  theory.  It  is  further  hoped 
that  the  present  formulation  can  be  used  as  a  stepping  stone  for  future  progress. 

2.2  Aerodynamics-Electromagnetics  Interaction  Equations 

The  rigorous  and  complex  governing  equations  for  flow  control  using  plasma  actuation  consist 
of  the  compressible  Navier-Stokes  equation,  Maxwell  equation  for  electromagnetics,  law  of  mass 
action  including  quantum  chemical-physics  for  nonequilibrium  chemical  reaction,  and  the  gas 
kinetics  for  transport  properties.  However,  from  all  existed  experimental  evidences  and  an  order- 
of-the  magnitude  analyses  [1-6,18-21],  the  electromagnetic  force  and  energy  appear  mostly  as  a 
perturbation  to  the  inertia  of  the  main  flowfield.  In  most  plasma  based  flow  control  environments, 
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the  Magnetic  Reynolds  number  which  is  the  product  of  the  electrical  conductivity  a,  magnetic 
penneability  pm,  characteristic  velocity  and  length  scales  of  the  study  phenomenon 

Rm  =  ul  /  ( o/lm )  ,  usually  has  a  value  much  less  than  unity.  Under  this  condition,  the  induced 

SB 

magnetic  flux  density,  —is  negligible  in  comparison  with  the  externally  applied  field  [4-6,19]. 

5t 

As  the  consequence,  the  Faraday’s  induction  law  can  be  decoupled  from  the  rest  of  the 
MaxwelFs  equations.  This  simplification  focuses  the  study  of  magneto-aerodynamic  interaction 
with  fluid  dynamics  rather  than  electromagnetic  wave  motion.  In  this  formulation  the 
electrostatic  force,  Lorentz  acceleration,  and  Joule  heating  are  appears  as  the  perturbing  source 
terms  in  the  generalized  Navier-Stokes  equations.  Therefore,  the  essential  physics  of 
aerodynamics-electromagnetics  interaction  for  plasma  flow  control  can  be  effectively 
approximated  by  a  simplified  governing  equation  system  which  can  be  summarized  as: 


^  +  V-[A(«  +  «,)] 
dt 


dwt 

dt 


(2-1) 


+  V  •  ( puu  +  pi  -  t )  =  peE  +  ( J  x  B ), 
dt 


(2-2) 


dpe 

dt 


+  V  ■  [peu  -  *57  T  +  Yj  PiU  A  +  qrad  +  u  •  pi  +  u  ■  f  ]  +  Qvt  -Qet=E-J. 


(2-3) 


The  vibrational  energy  conservation  equations  for  polyatomic  molecular  species  are; 

dpteiV  „  r  _  _  „  dw: 


-f-  +  V  •  [pt  ill  +  il  )eiV  +  qiV )]  =  eiV  ^  +  QV1 . 
dt  dt 


(2-4) 


The  electronic  energy  conservation  equation  has  been  traditionally  given  as; 

dp  e  -  dw ■  -  - 

—rf~  +  V  •  [p,.(t7  +  ui)ee  +  u-peI  +qe)\  =  ee—j-  +  E-J 
dt  dt 


+  [PeE  +  (d  x  B)]  ■  (u  +  uj )  +  Qe  s . 


(2-5) 

where  J  is  the  current  density,  J  -  e[-d+Vn+  +diVne  +d  V/?  +  (/?/./,  -nejue  -  n_jU_)E ]  of  the 
partially  ionized  gas,  and  can  be  approximated  by  the  transport  property  consisting  of  the  drift 
velocity  and  diffusion.  The  electron  and  ion  mobilities  ( ne,/u+ , //  )  are  traditionally  evaluated  as 
the  functions  of  the  reduced  electric  field  intensity  E/P,  and  the  diffusion  coefficients  {de,d+,d_) 

of  the  charged  molecules  can  be  calculated  by  the  Einstein  fonnula  [7,8].  The  diffusion  and  drift 
velocity  for  both  positively  and  negatively  charged  ions,  as  well  as,  electrons  are  included  to 
describe  the  individual  random  and  forced  motions  of  the  charged  particles.  The  vector  field  of  E 
and B  are  the  sum  of  the  externally  applied  and  induced  electrical  and  magnetic  field  intensities. 
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The  net  energy  transfer  between  translational,  vibration,  electronic  excitations,  and  the  sums  in 
each  degree  of  freedom  that  appears  in  Equations  (2-3),  (2-4),  and  (2-5)  are  designated  as  QVl,Qet, 

QV1 ,  and  Qe  l  respectively.  The  energy  cascading  processes  for  the  partially  ionized  gas  occur 

beneath  the  atomic  and  molecule  scales  which  cannot  be  described  the  gas  kinetic  theory,  thus 
must  be  modeled  through  the  quantum  chemical  and  chemical-physical  kinetics  [18-20]. 

It  is  important  to  note  that  in  the  governing  Equations  (2-2),  (2-3),  and  (2-5),  the  coupling 
between  fluid  motion  and  electromagnetic  force/energy  appears  as  the  electrostatic  force^f, 

Lorentz  acceleration  j  x  B  ,  and  Joule  heating  E-J  [6,7].  In  the  absence  of  an  externally  applied 
magnetic  field,  the  charged  particles  are  accelerated  by  the  electrostatic  force  due  to  charge 
separation  in  the  plasma  sheath.  The  well-known  induced  gas  motion  by  DBD  is  the 
consequence  of  momentum  transfer  by  collisions  between  charged  ions  and  the  neutral 
molecules  [1-3].  On  the  other  hand,  the  DCD  imparts  thennal  energy  to  the  external  stream 
through  Joule  heating  mostly  in  the  cathode  layer  in  addition  to  the  electrode  heating.  It  is  the 
basic  mechanism  of  plasma  actuator  by  thennal  effect  [3-5], 

The  vivid  illustration  by  applying  Lorentz  acceleration  to  alter  the  flowfield  structure  is 
displayed  in  the  classic  experimental  by  Ziemer  [22],  He  has  convincingly  demonstrated  that  an 
externally  applied  magnetic  field  can  change  the  bow  shock  standoff  distance  by  the  magnetic 
pressure,  b-b  /2ju  [19].  The  electromagnetic-aerodynamic  interaction  with  an  externally 
applied  magnetic  field  in  partially  ionized  gas  has  been  applied  to  hypersonic  control  surface  and 
inlet  [4,5],  a  sharp-nose  conical  configuration  by  Cristofolini  et  al.  [23],  and  heat-flux  mitigation 
at  the  stagnation  region  of  blunt  bodies  by  Matsuda  et  al.  [24]  and  Gulhan  et  al.  [25].  The 
dependence  of  the  interdisciplinary  fonnulations  on  the  Maxwell  equations  will  be  deferred  to 
the  section  of  the  discharge  modeling. 


The  law  of  mass  action  establishes  the  foundation  for  modeling  chemical  reaction,  and  the  rate  of 
species  production  and  depletion  in  the  species  conservation  law,  Equation  (2-1)  is  described  by 
the  Arrhenius  formula  through  the  forward  and  backward  chemical  reactions.  For  a  total  number 
of  the  chemical  species  of  n,  then  a  (n-1)  number  of  concentration  equations  are  required  to  be 
solved  simultaneously.  It  is  obvious;  the  required  effort  is  tremendous  and  tedious. 


dwt 

dt 


^2K->G 

7=1 


%  f  „  Yk'J 

II  {r 

k= 1  \Mk  ) 


Nlb  (  y*j 

-ui  5- 

4=1  \mk  y 


(2-6) 


For  most  investigations,  the  forward  and  backward  chemical  reaction  rates  are  adopted  from  a 
wide  range  of  chemical  kinetic  models  and  experimental  data  [13,15,19-21], 

The  definition  of  the  internal  energy  is  now  given  as; 


PC  =  Za(G/  +  +  YjPeeV,l  +ZM°  +Pe(Cv/e  +  (2-7) 


where  h’’  is  the  standard  heat  of  fonnation  for  all  reacting  species. 


To  be  consistent  with  the  kinetic  model  of  internal  structure  of  gas;  transport  properties  of  the 
gas  mixture  for  thermal  diffusion,  molecular  viscosity,  and  thenno  conductivity  need  to  be 
calculated  from  Boltzmann  equation  by  Chapman-Enskog  expansion  [26].  It  is  a  landmark 
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achievement  of  the  kinetic  theory  of  diluted  gas  mixture  by  its  ability  in  describing  the  transport 
property  of  any  combination  of  gaseous  mixture  by  the  inter-molecular  potential  functions.  The 
required  collision  integrals  and  cross  sections  have  been  obtained  by  either  the  Lenard- Jones 
potential  for  non-polar  molecules  or  a  polarizability  model  for  ion-neutral  non-resonant 
collisions  [27].  The  results  from  the  kinetic  theory  for  the  molecular  viscosity  and  thennal 
conductivity  of  individual  species  are  used  to  generate  a  global  property  for  a  gas  mixture. 

The  interdiscipline  computational  fluid  dynamics  (CFD)  equations,  Equation  (2-1)  through  (2-5), 
have  been  widely  applied  for  nonequilibrium  chemical  reacting  flows  and  are  independent  from 
the  specific  plasma  modeling.  Therefore,  the  numerical  algorithms  and  the  required  initial  values 
and  pertaining  boundary  conditions  will  not  be  elaborated  in  here.  The  details  can  be  easily 
found  in  open  literature  and  references  [4,5,13,15,18-21]. 

2.3  Physics-Base  Discharge  Modeling 

The  modeling  of  weakly  partial  ionized  gas  actually  consists  of  two  parts;  the  transport  property 
and  the  ionization  process  for  the  partially  ionized  gas.  The  physics  for  DCD  is  well  known 
through  the  luminary  contributions  by  Townsend  [7,8];  the  ionization  is  evolving  from  electrons 
impact,  secondary  emission,  and  avalanche.  The  physical  process  of  the  glow  discharge  has  long 
been  referred  to  as  the  Townsend  mechanism.  For  the  charged  carried  gas  mixture,  the  transport 
properties  must  include  the  force  diffusion  by  the  applied  electrical  field  on  an  electrically 
conducting  medium  and  the  ordinary  and  bipolar  diffusions.  For  the  gas  discharge  physics,  the 
modeling  and  formulation  based  the  drift-diffusion  approximation  for  the  transport  property  was 
established  in  the  late  1980s  [28,29].  The  model  for  DCD  flow  control  is  first  applied  by 
Surzhikov  and  Shang  and  extends  the  discharge  model  to  include  an  externally  applied 
transverse  magnetic  field  [30-32].  A  series  results  of  computational  simulations  using  this  model 
have  been  verified  and  validated  by  experimental  observations  conducted  in  a  Mach  five  plasma 
channel  [4,5,33]. 

For  supersonic  and  hypersonic  flow  controls  using  DCD,  both  electrodes  are  exposed  and 
flushed  mounted  over  the  surface  so  the  combined  electrode  and  joule  heating  will  raise  the 
temperature  of  the  gas  over  the  electrodes  locally.  In  turn,  the  elevated  local  temperature 
increases  the  displacement  thickness  of  the  shear  layer  to  induce  compressible  waves  and 
eventually  coalesces  into  an  oblique  shock  [4,5,30,31,34].  In  most  DBD  applications,  two 
electrodes  are  either  overlapped  or  placed  at  a  very  short  distance  of  a  few  mm  apart,  and  with  a 
minuscule  recessed  height  [1-3,10,1 1].  One  of  the  electrodes  is  exposed  and  the  other  is 
encapsulated  by  a  dielectric  coating.  The  most  significant  effect  of  an  encapsulated  DBD 
electrode  to  the  electromagnetic  field  is  by  modifying  the  electric  field  intensity  through  surface 
charge  accumulation.  Under  a  discharging  condition,  the  dielectric  electrode  becomes  a  capacitor 
to  diminish  the  electric  field  intensity  across  the  electrodes.  The  quenching  and  replenishing  rate 
of  the  surface  discharge  is  dependent  on  the  applied  alternating  voltage  across  the  electrodes  and 
the  surrounding  electrically  conducting  medium.  To  a  degree,  the  rate  is  also  dependent  on  the 
material  property  of  the  dielectric.  The  biased  electric  field  intensity  and  the  charge  separation 
over  the  electrodes  generate  a  periodic  electrostatic  force  to  move  the  charged  particles  which  in 
turn  transfer  momentum  by  colliding  with  the  neutral  particles.  In  fact,  the  momentum  transfer 
between  the  ions  and  the  neutral  particles  lead  the  so-called  electrical  wind  of  the  DBD.  It  is 
therefore  logic  to  model  the  partially  ionized  air  plasma  by  the  multiple-fluid  particle 
approximations. 
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In  DBD  operation,  the  Townsend’s  mechanism  still  controls  the  secondary  emission  initiated  by 
electron  collisions,  multiple  primary  avalanches,  and  ultimately  maintains  the  discharge. 
Impressive  pioneering  efforts  by  Massines  et  al.  [35]  and  Lee  et  al.  [36]  have  gained  helpful 
insights.  All  the  aforementioned  approaches  are  assumed  that  in  a  high-pressure  discharge  field 
the  electrons  and  ions  are  in  equilibrium  with  the  electric  field.  In  particular,  Lee  et  al  solve  the 
Boltzmann  equation  for  the  electron  energy  distribution  function.  All  previous  investigators 
understood  the  key  elements  of  the  discharge;  the  temporal-spatial  variation  of  the  physics  must 
be  analyzed  simultaneously.  Both  efforts  investigated  the  basic  phenomenon  during  a  single  AC 
cycle  after  the  initial  transient  has  subsided. 

The  individual  velocities  of  charged  particle  motions  of  weakly  ionized  plasma  that  constitutes 
of  neutral  molecules  (10  cm'  ),  positively  (10  cm"  )  and  negatively  charged  ions  (10  cm"  )  , 
as  well  as  electrons  have  been  derived  by  Surzhikov  et  al  [3 1,32]  from  the  momentum  equation 
of  charged  particles  motion.  The  electron-neutral  collision  frequency  is  proportional  to  the 
charged  particle  number  density  and  the  electron  temperature.  The  time  scale  is  orders  of 
magnitude  shorter  than  the  organized  motion  of  the  neutral  particles;  thus  the  time  dependency  of 
the  molecular  shear  stress  and  inertia  of  neutral  particles  in  the  momentum  equation  of  charged 
particles  motion  are  negligible.  By  defining  the  partial  pressure  of  each  charged  species  and 
omitting  the  relatively  slower  velocity  of  the  organized  motion  of  neutral  particles,  the  individual 
charged  species  velocities  can  be  summarized  as  the  following; 

neue  =  -deVne  -  / ueneE  -  /uene(ue  x  B) 


n+u+  =  -d+Vn+  +  /J+n+E  -  ju+n+(u+  x  B)  (2-8) 

n_u_  =  -dVn_  -  / u_n_E  -  ju_n_(U_  xB) 


It  is  important  to  recognize,  the  above  equations  actually  define  the  diffusion  velocities  of  an 
inhomogeneous,  electrically  conducting  mixture.  The  first  terms  in  Equation  (2-8)  describe  the 
ordinary  diffusion  due  to  local  concentration  of  electrons,  positively  and  negatively  charged  ions. 
The  second  terms  are  the  force  diffusion  in  form  of  drift  motion  for  an  electrically  charged 
species  by  an  externally  applied  electric  field.  The  last  tenns  are  also  the  force  diffusion  by  the 
Lorentz  acceleration. 


The  species  conservation  equation  becomes  a  simplified  approximation  to  the  continuity 
equation,  Equation  (2-1)  for  plasma  flow  control.  The  diffusion-drift  model  of  weakly  ionized  air 
plasma  is  equally  applicable  to  DCD  and  DBD.  This  formulation  has  been  adopted  by  most 
recent  numerical  simulations  by  Surzhikov  et  al  [30-32],  Gibalov  et  al.  [37],  Boeuf  et  al.  [38-39], 
Unfer  et  al.  [40],  Solv’ev  et  al.  [41],  Likhanskii  et  al.  [42],  Shang  et  al  [43-45],  and  Huang  et  al. 
[46].  One  of  the  weaknesses  of  this  approximation  is  that  the  detailed  distinction  between 
metastable  molecules  is  no  longer  recoverable  by  this  approximation  [9-11].  The  weakly  ionized 
gas  is  treated  as  a  medium  consists  of  three  charged  and  one  neutral  species.  The  species 
conservation  equations  for  the  three  distinct  charged  particles  are: 


dt 


+  v-fe 


dwe 

dt 


(2-9) 
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(2-10) 


dn+ 

dt 


+  vr 


dw+ 

dt 


dn 

dt 


+  V-f 


dw 

dt 


(2-11) 


dwe  dw+  dw_ 

Where  the  terms  in  the  right-hand-side  of  the  equations;  — —  ,  — , —  ,  and  — , —  designate  the 

at  dt  dt 

production  and  depletion  rates  of  the  considered  ionized  species.  These  rates  have  either  been 
determined  by  some  chemical  kinetic  models  [9,13,15,41,47]  or  simplified  approximations  based 

on  the  known  ionization  process  [30-32,38-46],  The  symbols  Te,r+,  and  T  denote  the  electron, 
positively  and  negatively  charged  ions  flux  densities  and  have  the  following  definitions; 

? e  =  -  dyne  -  ne/ae{E  +  uexB ) 

=-d+Vn+  +n+Ju+(E  +  u+  xB)  (2-12) 

=  -dyn_  -  n_ju_(E  +  u_x  B ). 


The  important  connection  between  the  coefficients  of  diffusion,  dt  and  the  mobility,  //,  of  drift 
velocity  is  given  by  the  Einstein  relation  [7,32],  the  relative  magnitude  is  a  function  of  the 
electron  characteristic  energy  to  the  elementary  electric  charge. 


kT 

Be  =  ~  Be(  E  /p);  de  =  — *-/ie 

meven  '  e 

€  kT 

H±  /A  (  E  / P^) ,  d±  ~  B± 

m±v±n  1  e 


(2-13) 


Like  the  energy  spectrum,  it  is  a  function  of  the  reduced  electrical  field  |j?|  /P.  When  a 

considerable  free  space  charge  is  formed  by  charged  separation,  the  polarized  electric  field 
restrains  the  motion  of  electrons  from  that  of  ions  to  appear  as  the  ambipolar  diffusion.  This 
coefficient  of  diffusion  has  frequently  approximated  by  a  function  of  the  ratio  of  the  electron 
energy  and  ambient  states  [30-32], 


A  general  formulation  for  the  drift-diffusion  model  including  an  externally  applied  magnetic 
field  is  very  complex  because  the  gyration  of  charged  particle  associated  with  the  magnetic  field 
polarity  [6].  For  the  purpose  of  an  analytic  formulation,  it  is  convenient  to  resolve  the  velocity  of 
charged  particles  into  component  parallel  and  perpendicular  to  the  direction  of  an  externally 
applied  magnetic  field.  In  practical  application,  the  transversal  magnetic  field  for  flow  control  is 
known  to  be  the  most  effective  [6-8].  Therefore,  a  two-dimensional  formulation  in  (x,y)  plane  for 
charged  particle  flux  densities  including  all  the  electromagnetic  effect  is  viable  by  restricting  the 

magnetic  flux  density  to  the  z  coordinate  B , .  The  resultant  equations  are  similar  to  that  in  the 


absence  of  an  externally  applied  magnetic  field  by  introducing  effective  electric  field  strengths  as, 
b  /•;  /•;.  bE.  +  Ev 

— ,  and  the  components  of  the  electron  flux  density  becomes; 


Ee,x  = 


1  +  bt 


’Ee,y=- 


1  +  bt 
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r  -_„nv _ 1  A  dn*  i  b'  A  dn~ 

1  e,x  Mene^e,x  i  ,  i2  ®e  -  ,  d„ 


1  +  bze  e  dx  1  +  Zr  P  Sy 


r  =  -„«£  — ‘ 

^  Me  *  ^  ]+h2  e  a..  ,  ,  L2  « 

e 


\  +  b 2  e  dy  \  +  b 2  e  8x 


Ex  +  b+E 

Similarly  for  the  positively  charged  ions,  by  letting  E+x  = - ~r—,E 


1  +  b 


2  ’J 


r  _  „  „  F _ Lj  _ a" 

*  +  V  MMK8+.x  ,  .  ,  2  %  -  ”  ^ 


1  +  b2  Sx  1  +  b2  +  <3y 


1  ,  <3/2  h  ,  <3/2 


r ,  v  =  ju.n.E ,  v - 7<i+  — -  +  — ^c/ 

+,>*  '  +  +  +,r  1  .  r  2  +  ^  1  .  r  2 


1  +  b2  dy  1  +  b2  +  Sr 


27  = 

-■*  1  .  7.2  ’  -J- 


Ex-bEv 

Finally,  set  the  effective  electric  field  strength  as  E  x  -  — - — — 

components  of  the  negatively  charged  ions  are; 

j-h  r  \jdn_  b_  dn_ 

F  =  -u  n  E - Td  - + - -d  - 

— 1  i  +  h:  '  &  i  +  b:  -  dy 


r  \  dn_  b_  dn_ 

T  v  =  -H_n_E_^  Yd_  - - —d_ 


\  +  b_  dy  1  +  b_  dx 


(2- 14a) 

(2- 14b) 


+,y 


Ev  ~  b+E  x 
1  +  b; 


,  we  have 


(2- 15a) 

(2- 15b) 


i  +  h2 


.  The 


(2- 16a) 

(2- 16b) 


In  the  above  equations,  b  ,b+,  and  /,  are  the  Hall  parameters  for  electrons,  positively  and 
negatively  charged  ions,  which  are  directly  related  to  the  Lannor  frequencies  of  electrons  and  the 


ions, 


b  = 


Me8, 


=  oy  b  =  m_Bz 


u 


c 


=  and  b+  =  2^— i  * 
V  C  V 


—  •  In  here,  ve ,  v_  ,  and  v  are  the 


averaged  electron  and  ion  collision  frequencies  [6,7]. 


The  above  fonnulations  have  been  validated  by  successfully  simulating  hypersonic  flow  controls 
in  the  presence  of  an  externally  applied  transverse  magnetic  field  at  a  sharp  leading  edge  and  in  a 
virtual  variable  area  inlet  cowl  [4,5,30,33],  The  numerical  results  are  verified  by  experimental 
data  generated  by  a  Mach  five  plasma  channel.  On  the  sharp-leading  edge,  the  experimental  data 
and  numerical  results  of  the  surface  pressure  reveal  a  substantially  amplified  viscous-inviscid 
phenomenon  by  the  transverse  magnetic  field  through  the  magneto-aerodynamics  interaction. 
Similarly,  the  computational  simulations  also  compare  well  with  the  Pitot  pressure  surveys  in  the 
cross-section  planes  and  surface  pressure  distributions  within  a  square  cross-section  inlet  [4,5]. 


2.4  Ionization  Modeling 

The  ionization  process  by  electron  impact  is  distinct  from  the  thennal  excitation,  but  the 
nonequilibrium  ionized  species  concentrations  must  satisfy  the  law  of  mass  action.  For  this 
reason,  the  ionization  processes  are  treated  as  chemical  reactions  of  molecules,  atoms,  and 
electrons  or  ions.  Chemical  kinetics  models  for  the  low-temperature  plasma  generation  span  a 
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very  wide  range  of  complexity,  from  the  detailed  description  by  Elisson  et  al.  [9],  Bogdanoff  et 
al  [13]  to  simplified  approximations  by  Solov’ev  et  al.  [41]  and  Singn  and  Roy  [47].  Singn  et  a.l 
reduce  the  chemical  kinetic  model  by  omitting  the  metastable  species  and  the  nitrous  oxide. 
Solov’ev  et  al.  simplify  their  selection  of  the  elementary  chemical  reactions  based  on  the 
comparable  reaction  rates  of  species  to  the  characteristic  time  scale  on  the  order  of  30 
nanoseconds  for  DBD  in  atmospheric  air.  The  approach  is  rational  and  the  fidelity  of  numerical 
simulation  to  nonequilibrium  chemical  reaction  is  solely  depended  on  the  empirical  detennined 
rate  constants.  Unfortunately,  it’s  the  weakest  link  of  ionization  models  for  both  the  high- 
temperature  thennal  excitation  and  electron  collision  mechanisms  [9,13,20,21,41],  The  complete 
differential  equations  system  of  species  conservation  is  very  stiff  to  become  a  fonnidable 
computational  challenge. 

An  alternative  approach  in  determining  the  net  rates  of  generation  and  depletion  of  individual 
ionized  species  can  be  achieved  by  the  chemical-physics  ionization  modeling.  In  fact  the 
ionization  by  electron  collisions  occurs  at  the  outer  limited  of  the  Maxwell-Boltzmann 
distribution  where  the  gas  energy  is  lower  than  the  ionization  potential  [7],  The  dominated 

species  of  DBD  are  the  metastable  molecules  N2  and  C>2(  bl1?d  ).  The  main  ionized 

processes  of  a  nonequilibrium  volumetric  discharge  consist  of  four  type  of  reactions; 
electron/molecular,  atomic/molecular,  decomposition  and  synthesis  [7,9].  At  any  instance  during 
discharge,  the  ionized  species  concentration  is  the  net  balance  between  ionization,  detachment, 
attachment,  and  recombination  processes.  The  ionization  model  is  therefore  focused  on  these 
mechanisms. 

The  classic  formulation  for  electronic  impact  ionization  is  the  similar  law  by  Townsend  and  is  an 
empirical  formula.  In  the  formulation,  the  ionization  coefficient  «  which  measures  the  number 
of  ionization  by  electron  per  unit  distance  is  a  function  of  the  reduced  electrical  field  |j?|/P  [7,8]. 

This  quotient  is  also  a  measure  of  the  energy  gain  by  a  charged  particle  between  collisions  from 
the  principle  of  similarity.  In  fact,  the  coefficient  a  =  A  exp  (-Bp  /  |a|)  of  Townsend’s  similarity 

law  holds  extremely  well  in  comparison  with  a  large  group  of  experimental  data  both  in  the 
ionization  frequency  and  degree  of  ionization.  For  discharge  in  air;  A=15  and  B=365  in  the  E/P 
range  from  100  to  800  (V/Torr  cm).  At  a  relatively  higher  value  of  E/P,  an  accuracy 
improvement  may  be  needed  but  is  not  essential.  The  coefficients  of  ionization  for  air  can  be 
summarized  as; 

(— )  =  15exp[ — r/77— — ],  1/  cmtorr  (2-17) 

P  (| E\/P ) 

The  dissociative  recombination  is  the  fastest  mechanism  of  the  bulk  recombination  of  the  weakly 

ionized  gas,  and  is  a  simple  binary  chemical  reaction.  The  decay  rate  with  time  of  plasma  is  often 

dfi  7  3 

given  as  ( — ^).  =  -J3n+ne ,  and  typically  coefficient  (3  is  assigned  a  value  of  2x10'  cm  /s  and  the 
dt 

3 

characteristic  decay  time  scale  is  less  that  10'  second. 

Another  main  mechanism  of  charge  neutralization  is  the  ion-ion  recombination  process.  At  the 
low  pressure  environment,  the  recombination  takes  place  through  binary  collision  and  the 
reaction  is  similar  to  charge  transfer  [9,13,15].  At  the  moderate  pressure,  the  process  proceeds 
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through  triple  collisions,  again  the  rate  of  ion-ion  recombination  can  be  given  as  =  ~P,n+n_  ■ 

dt 

For  an  example,  the  recombination  rate  constants  between  02  andff ,  as  well  as,  NO  and  N02 

25  26  6 

are  on  the  order  of  magnitude  of  10'“  and  10'“  cm  /s  [7,13]  .  The  maximum  value  of  the  ion-ion 
recombination  at  the  one  atmosphere  has  a  value  of  10'  cm  /s  [7]. 

The  electron  attachment  and  detachment  are  the  formation  and  depletion  of  the  negative  charged 
ions  in  the  partially  ionized  air  [9].  For  the  strongly  bound  molecules,  a  sufficiently  high  energy 
is  required  for  the  dissociative  attachment  of  an  electron.  The  electron  attachment  is  the  main 
mechanism  for  removing  the  electron  from  the  negatively  charge  ions.  The  lost  electron  number 

density  can  be  given  as;  (— )a  =  -v  n  ■  The  attachment  frequency  of  electron  in  dry  air  at  one 

dt 

o 

atmosphere  condition  is  around  10  /s  [7]. 

The  coefficients  of  the  detachment  of  electron,  Kd  in  partially  ionized  gas  at  the  room 
temperature  have  the  value  of  10'10  cm3/s.  The  value  for  the  metastable  molecules  Ni(  )  and 

02(b%)  in  air  are  unknown  [9,13].  But  by  an  indirect  estimate,  the  discharges  are 

characterized  by  a  value  of  10‘14  cm  ’/s.  It’s  also  interesting  to  note  that  the  weakly  ionized 
plasma  is  sustainable  in  negatively  charged  gas  at  a  lower  value  of  E/p  than  that  by  a  short  pulse 
discharge. 

From  the  above  brief  discussion,  the  plasma  generation  and  depletion  processes  are  modeled  as 
the  following  [30-32]; 


bne 

dt 

dn+ 

dt 

dn_ 

dt 


~  V  •  (ne/ueE  +  deVne )  =  a(\E\)  \  fe  \ -f3n+ne  - v ane  +  Kdnnn_ 


+  V  •  ( n+ju+E  -  d+Vn+ )  =  a(  E  )  |  f  |  -j3n+ne  -  /3in+n_ 


-V  •  (n_/j_E  +  d_Vn_)  -  v  ane  -  rcdnnn_  —J3in+n_ 


(2-18) 


In  fact,  the  sum  of  the  above  equations  satisfies  fully  the  charge  conservation  law  which  is 
derived  from  the  Maxwell  equations  [50].  Again  in  the  above  approximation;  the  number  density 
of  neutral  particles  of  air,  «  is  treated  as  a  constant  to  have  a  value  of  2.69x  1 0  /cm  .  The 

attachment  frequency  vy  is  formally  defined  as  v  =  («  /  p)uedrift  ■  p  ■  In  here,  a  commonly 
adopted  value  is  (aa  /  p)  =  0.005  /  cm  ■  ton ■ ,  and  ue  dn,t  is  the  drift  velocity  of  electron.  From  the 
drift-diffusion  theory,  the  velocity  can  be  consistently  evaluated  as  ue  drjft  =  //  //I . 


The  coefficient  of  ion-ion  recombination  //  is  estimated  from  the  chemical  kinetics  of  ionized 
and  molecular  nitrogen  and  oxygen  to  have  a  value  of p.  =  i.6  x  10  7 , end  /  s  •  The  detachment 
coefficient  is  estimated  to  have  a  range  of  io  14  <  Kd  <  8.6  x  io10,cm3  /  5-  •  This  coefficient  is 

important  for  the  negatively  charged  ion  number  density  computation;  because  it’s  the  principal 
mechanism  that  governing  the  depletion  of  this  species.  Fortunately,  its  value  can  be  verified  by 
the  discharged  electrical  current  from  experimental  data. 
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The  above  model  of  electron  collision  ionization  is  actually  replaced  the  energy  conservation 
equations  of  the  vibrational  and  electronic  internal  excitations,  Equations  (2-4)  and  (2-5),  as  well 
as,  the  law  of  mass  action  Equation  (2-6)  for  chemical  kinetics  modeling.  The  electron 
attachment  mechanism  has  been  modeled  to  be  proportional  to  the  Townsend’s  ionization 
process  for  a  non-physical  but  computational  stability  concern.  In  this  approach,  a  portion  of  the 
negative-ion  ionization  is  split  from  the  electron  generation  [40], 

In  theory,  the  above  ionizing  models  can  be  further  improved  by  expanded  the  experimental  data 
base  or  by  ab-initio  approach  to  the  quantum  chemical-physics.  In  this  regard,  the  state-of-the-art 
status  in  low-temperature  ionization  modeling  is  similar  to  that  of  chemical-physics  models  by 
the  approximated  chemical  reaction  rates  [13,41,47],  but  it  has  a  decided  advantage  by 
substantially  reducing  computational  resource  required  for  simulations. 

2.5  Electromagnetic  Equations 

The  Maxwell’s  equations  govern  all  electromagnetic  phenomena  and  consist  of  the  Faraday’s 
induction  law,  the  generalized  Ampere’s  circuit  law,  and  two  Gauss’  divergence  equations  for 
magnetic  flux  density  as  well  as  electric  displacement  [50].  For  the  DCD  and  DBD  modeling,  the 
governing  equations  including  the  boundary  conditions  are  based  on  these  fundamental  laws  of 
electromagnetics.  The  differential  equations  system  is; 


SB 

dt 

dD 

dt 


+  Vx£  =  0 
—  VxH  =  —J 


VB  =  0 


(2-19) 


V  ■  D  =  p 


The  associated  boundary  conditions  of  the  differential  system  bounded  by  media  interface  can  be 
summarized  as; 


n  x  (Et  —  E2  )  =  0 

n  (A  -A)  =  Ps 

n  x(Hl-H2)  =  Js 
n  ■  (Bx  —  B2  )  =  0 


(2-20) 


The  above  boundary  conditions  require  the  parallel  component  of  electric  field  strength  E  on  a 
media  interface  to  be  unaltered,  but  the  difference  in  the  electric  displacement  component  D 
across  the  interface  must  be  balanced  by  the  surface  charge  ps.  At  the  same  time,  the  magnetic 
field  strength  H  parallel  to  the  medium  interface  needs  to  be  balanced  by  the  electric  surface 
current  Js,  and  the  normal  component  of  the  magnetic  flux  density  B  to  the  media  interface  is 
unchanged. 

In  most  plasma  based  flow  control  environments,  the  Magnetic  Reynolds  number  usually  has  a 
value  much  less  than  unity  [4,6,19].  In  the  absence  of  an  externally  applied  magnetic  field  its 
induced  magnetic  flux  density  is  negligible  and  the  electric  field  dominates  the  aerodynamic 
interactions  of  DCD  and  DBD.  This  simplification  reduces  the  electromagnetic-aerodynamic 
interaction  equations  into  the  weak  form.  In  this  formulation  the  electrostatic  force,  Lorentz 
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acceleration,  and  Joule  heating  are  appears  as  the  source  terms  in  the  interdisciplinary  equations, 
Equation  (2-2)  and  (2-3). 

In  deriving  the  species  concentration  equation,  the  Maxwell  equation  is  not  explicitly  include  but 
the  sum  of  Equations  (2-2),  (2-3),  and  (2-4)  fully  satisfies  the  charge  conservation  equation.  This 
equation  is  obtained  by  taking  the  divergence  of  the  generalized  Ampere’s  circuit  law  and 
invoking  the  Gauss’s  law  for  the  divergent  free  of  the  magnetic  flux  density  [50]; 

%  +  V-J  =  0  (2-21) 

ot 


In  the  above  equation,  the  electric  charge  density  and  the  electric  current  density  are  defined  as 
[30-32]; 


pe=e(n+-n_-ne) 

J  =  c(f+-f_-fe) 


(2-22) 


The  DCD  and  DBD  plasma  model  is  formally  closed  by  the  Poison  equation  of  plasmadynamics 
which  satisfies  the  Gauss’s  law  for  the  electric  displacement,  v  •  (eE)  =  p  and  assumes  the 
electric  field  intensity  can  be  derived  from  the  electrical  potential.  Specifically  the  electrical  field 
intensity  equals  to  the  negative  gradient  of  the  electric  potential,  e  =  -\/<p  .  The  resultant  well- 
known  Poisson  equation  appears  as: 


V  V  =  —  (ne  +  n_  -  n+ )  (2-23) 

s 

The  appropriate  boundary  conditions  for  DCD  and  DBD  simulations  associate  with  the  unique 
characteristic  of  different  discharges  will  be  discussed  in  the  following  sections. 

It  is  recognized  that  the  DBD  modeling  may  not  necessarily  be  adequate  to  describe  the  plasma 
sheath  with  free-space  discharge  separation,  especially  in  conjunction  with  the  Poisson  equation 
of  plasmadynamics.  The  shortcoming  is  the  consequence  that  the  compatible  electric  field  is 
evaluated  from  the  basic  assumption,  in  that  the  electric  field  intensity  is  directly  derivable  from 
the  electric  potential.  Numerous  attempts  have  been  made  to  improve  this  deficiency  of 
plasmadynamics  equation  [26,27,48]  including  a  Monte-Carlo  simulation  by  Boeuf  et  al  [49]. 
Although  these  improvements  have  captured  some  new  details  of  the  gas  discharge  phenomenon, 
but  the  basic  approximation  is  still  capable  to  predict  all  essential  features  with  sufficient 
accuracy  for  the  motion  of  partially  ionized  particles. 

2.6  Numerical  Procedure 

The  governing  equations  for  DCD  and  DBD  modeling  are  identical  and  the  distinctive  physics 
are  generated  by  the  different  initial  values  and  boundary  conditions.  The  physics-based  model 
consists  of  a  mixed  inhomogeneous  hyperbolic  and  elliptic  partial  differential  equations, 
Equations  (2-18)  and  (2-23).  The  equations  are  source  terms  dominant  and  the  eigenvalues  of 
this  differential  equations  system  span  a  large  range  to  make  it  very  stiff  for  computation.  A 
number  of  innovative  measures  must  be  introduced  to  ensure  stability  and  accuracy  of  the 
numerical  simulation.  In  the  earlier  efforts,  the  entire  equations  system  was  solved  by  a 
straightforward  implicit  scheme  which  possesses  a  narrow  convergent  band.  Therefore  a  set  of 
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stringent  initial  values  must  be  specified  to  meet  a  converged  tolerance  within  a  reasonable  time 
[31,32]. 

The  most  recent  approach  implements  the  total  variation  diminishing  (TVD)  discretization  by  an 
implicit  time  accurate  marching  scheme  coupled  with  multigrid  method,  and  the  implicit 
boundary  conditions  on  electrodes.  The  numerical  procedure  attempts  to  establish  a 
computational  stable  and  efficiency  method  for  simulating  surface  DBD  plasma  actuator.  First, 
the  two-dimensional  charged  species  conservation  law,  Eq.  (2-18),  is  discretized  over  the  control 
volume  as  shown  in  the  following  for  the  variable  n+ 


dn+  _  _^(n+E x\i+l/2,j)  (n+^x\i-l/2,j)  +  (W+^y)(iJ+l/2)  (n+^y  )(i  J-l/2)  - 


dt 


1  r/  ,  dn  dn+ 

+  dx  (,+1/2J)  ^ 


4 


0 i,j ) 


(2-24) 


Ar 


(i,j) 


1  dn+ 

+ - [«— 


V,,,,''  dy'VJ+U2)  '  '  dy 

In  the  above  formulation.  The  diffusive  tenns  were  treated  with  the  central  differencing  and  the 
source  terms  are  separated  in  two  parts;  for  the  term  s,.  does  not  explicitly  contain  the  pertaining 

dependent  variable  n  and  the  nonlinear  termS  is  a  product  of  collision  process  with  other  charged 
species.  The  typical  convective  flux  on  the  control  surface  («+je;v)(.+1/2  can  be  expressed  as: 


dx 

«— )(i>/— 1/2  )]  +  s+-sn+ 


(n+^x\i+l/2,j)  (Ex)(m/2 jj  [(”+)(/ J)  +(n+)(WJ)]/2+  I  {^x)(i+y2 J)  I  KA”+)(/,/)  +(A”+)(/+lJ)]/4 

“  (^)(,+l/2  J)  +  )(i+V2  J)  [(A”+)(7+l,y)  _(A”+ )(;,;)]/[(”  +  )(<+lJ)  -K)(u)]| 

[(«+W-,-K)(/j)]/2 


The  TVD  limiter  by  Yee’s  MINMOD  scheme  [51]  is  adopted  as; 


(2-25) 


M(U)  =  V/>  max[0,min(|  A n(MJ)  |,s((J)An(,.u)]  (2-26) 

where  /\/7  =  w  —  yi  and  ?  =  si qvl  ( Aw  ) . 

The  discretized  charged  species  conservation  equations  are  expressible  by  a  simple  algebraic 
equation  of  the  following  form; 

dn+,UJ)  _  „  „ 

Qt  aE,(i,j)n+(.i+hj)  +  +  aN,(i ,j+l)n+OJ+l)  +  aS,(iJ)n+(iJ- 0  (2-27) 

-(aPXU)+sV.J)>HU)+s+.«J) 

One  more  step  we  took  in  this  approach  is  to  use  the  Delta  formulation  via  the  iteratively 
diminishing  residual  approach.  By  setting gn+  =  n"  - n\ ,  where  **  is  the  current  iteration  level 
and  *  is  the  previous  iteration  level,  the  discretized  equations  system  can  be  summarized  as; 
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(2-28) 


8Sn+(i  n 

=  aE,(i,j)^nHM,j)  +  aw,(i,j)^n+(i-i,j)  +  aN,(i,j+D^n+(i,j+\)  +  as,(i,j)^nm,j- 1) 
~(ap,u.i)  +  s(,,j))dn+(u)  +  RHS+(i  j) 


In  Equation  (2-28),  the  RHS+  n  is  an  abbreviation  of  the  right-hand-side  of  Eq.  (2-24),  or  Eq. 

(2-27),  according  to  the  diminishing  residual  delta  formulation  to  be  evaluated  at  the  previous 
iterative  level  *. 

The  third-order  TVD  type  of  Runge-Kutta  algorithm  for  time  advancement  is  applied  by  the 
consecutive  implicit  operators.  The  following  formulations  contain  the  RHS  of  Equation  (2-27) 
at  different  fractional  time  steps  and  the  previous  time  level  to  appear  as; 


1 


[^t+ap,(uj)+S(ij)]8n+,(ij)  ~  aE,(ij)^n+,mj) +  awxij)^n+,a-ij) + 


aN,(i,j)Sn%,j+l)  +  aSJLt,J)Sn%,J- 1)  +  RHS-U,n 


Yl  —  Tl*1  +  SyI^ 
n(iJ)  ^  UndJ) 


8n\)\,  + 


[—  +  +  Sfn  }8n%)  =  -  -  „„ 

At  4  °’J)  U’J)  4  At  (,J) 

-^(aE,(I,j)dn+Xi+l,j)  +  aw,(i,j)^n+ll-hj)  +  aN,(i,j)^n+li,j+ 1)  +  as,(i,j)dn+,(i,j-l)^  + 

nlj)=nlj)+Snlj) 

[—  +  ~aP(ii)  +  S <2>  ]Sn "+‘  =  — —  Sn™  + 

At  3  X’J)  (J)  (J)  3  At  (’J) 

^(aE,(IJ)8KxMJ)  +  aw,(iJ)SnT,(I-\,j)  +  aN,(i,j)^n+XiJ+ 1)  +  aS.(i,j^n'Z),j- 1))  +  3 ~RHS+l,J) 

,„«+l  ,2  .  c*  n+ 1 

Yl,  x  —  Yl,  x  +  OYl,  x 

U,j)  0,j )  U,j) 


(2-28) 


Similarly,  the  Poisson  equation  of  plasmadynamics,  Eq.  (2-23),  can  be  discretized  using  the 
central  differencing  scheme.  In  the  Delta  formulation  the  equation  yields; 


lp,l 


u  Mi 


j) 


a 


8(p{i 


'(/+!  J)  +  aW,(‘J) 


( i-ij )  +  aN,(ijy 


8cp. 


\U+ 1)  +  aSXi,jMiJ-l)  + 


(2-29) 

where  RES  is  the  residual  value  of  Eq.  (2-23),  namely  the  RHS  subtracting  LHS.  Both  Equations 
(2-28)  and  (2-30)  require  a  very  efficient  implicit  solver.  The  SIP  (Strong  Implicit  Procedure)  or 
the  ILU  (Implicit  Lower  and  Upper)  decomposition  scheme  by  Stone,  or  a  modified  version  of  it, 
MSIP,  is  adopted  [52,  53].  The  algorithm  is  very  efficient  for  solving  a  sparse  matrix  equation 
system;  AU  =  R  .  In  here,  U  may  be  charge  species;  ne , n+ ,  n_  or  electric  potential,  (p;  A  is  the 

five-point  stencil  matrix  whose  coefficients  are  made  up  from  a '  s  shown  in  Eq.  (2-29)  and  (2- 
30).  By  addition  a  complementary  matrix  C  to  the  A,  the  resulting  matrix,  M,  can  be  solved 
implicitly  using  the  lower  and  upper  matrix  decomposition: 
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M  =A  +  C 


(2-30) 


The  identical  implicit,  iterative  solving  procedure  is  then; 

MU**  =  CU*  +  R  (2-31) 

Stone  further  introduced  a  partial  cancellation  term  in  both  side  of  the  Eq.  (2-32)  in  order  speed 
up  the  rate  of  convergence.  The  readers  are  referred  to  the  original  paper  for  the  details  of  the 
procedure  [52,53]. 

For  the  surface  discharge  phenomenon,  the  steep  gradient  always  occurs  near  the  comer  of  the 
exposed  electrodes  and  dielectric  surface,  thus  a  high  non-uniform  mesh  is  generally 
recommended  to  solve  the  equations.  Due  to  a  large  number  of  highly  stretching  grid  points 
were  used  to  cover  the  whole  computational  domain,  the  multigrid  method  offers  the  best  solving 
option  [54].  The  multigrid  method  renders  the  error  of  the  rapidly  varying  Fourier  component  to 
be  first  eliminated  by  the  fine  mesh  system.  The  slowly  decay  error  associated  with  the  spectral 
radius  becomes  a  smooth  function.  Through  this  implementation,  an  amazing  convergent 
acceleration  is  observed,  but  the  most  important  improvement  to  the  numerical  simulation  is  the 
computational  stability  by  the  combination  of  SIP  (or  MSIP),  the  implicit  third-order  Runge- 
Kutta  scheme  and  TVD  spatial  discretization,  which  contributes  greatly  to  the  accurate  of  the 
numerical  solution. 

In  Figure  1,  the  comparative  study  of  the  implicit  solving  schemes  is  depicted  on  a  highly 
stretched  (641  by  891)  grid.  The  remarkable  convergent  acceleration  from  the  conventional 
tridiagonal  line  relaxation  (TDMA)  to  SIP  and  finally  the  combination  of  SIP  and  multigrid 
procedure  is  clearly  illustrated.  The  acceleration  in  residual  reduction  by  the  SIP  algorithm  alone 
exceeds  four  orders  of  magnitude.  The  rate  of  convergence  is  further  vastly  improved  by  the 
combined  application  of  the  SIP  and  multigrid  techniques,  although  the  computational  stability  is 
not  easily  quantified  but  its  impact  is  impressive.  Through  these  superior  computational 
attributes,  the  grid  independent  solutions  are  truly  realizable. 
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Figure  1  Iterative  Convergent  Acceleration  for  Solving  Poison  Equation  of  Plasmadynamics 

Figure  2  demonstrates  the  ability  in  achieving  the  grid  independent  solutions  through  the 
combination  of  SIP  and  multigrid  technique.  The  electric  field  potential  vectors  are  projected  on 
the  contours  of  the  charge  density  distribution.  A  total  three  different  nonuniform  grids  of 
(257x417),  (641x897),  and  (1281x2049)  are  used.  In  the  left  side  of  the  graph,  the  exposed 
electrode  is  positively  biased,  and  the  counter  part  of  the  negative  polarity  is  at  right  side.  All  the 
solutions  are  interpolated  to  a  fixed  position  for  purpose  of  comparison.  After  the  second  grid 
density  refinement,  the  difference  between  solutions  is  undiscemible. 
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Figure  2  Grid  Independent  Studies 

(a-1)  to  (a-3)  are  for  positive  cycle  and  (b-1)  to  (b-3)  are  for  negative  cycle.  Contours  show  the  charged 
density  distribution  and  vectors  are  the  force  field.  The  solutions  were  interpolated  to  fixed  position  for 

purposes  of  comparison. 


2.7  Boundary  Conditions  of  Direct  Current  Discharge 

The  effect  of  DCD  in  applications  to  flow  control  is  derived  from  the  dissipative  Joule  heating  of 
the  discharge  [4-6].  However,  the  thermal  effect  also  includes  the  electrode  heating  from  the 
metallic  electrodes  and  mostly  over  the  cathode.  The  electrodes  generally  have  two  typical 
configurations;  the  parallel  electrodes  separated  by  a  gap  distance  d,  and  the  side-by-side 
arrangement.  The  boundary  condition  for  DBD  is  straightforward: 

On  the  cathode,  the  electric  field  potential  vanishes  on  metallic  surface,  all  negatively  charged 
ions  are  repelled,  and  the  number  density  of  the  positively  charged  ions  reaches  the  cathode 
unaltered.  The  secondary  electron  emission  is  the  result  of  the  Penning  tunneling  through  a  very 
strong  electrical  field  intensity  in  the  cathode  layer  generated  by  the  positive-ion  accumulation 
immediately  adjacent  to  the  solid  surface.  The  quantification  of  the  electron  emission  from  the 
cathode  by  positively  charged  ions  is  described  by  an  emission  coefficient  y  which  has  a  well- 
known  value  range  from  0.01  to  0.1,  depending  on  the  electrode  material  [7,28,29] 
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(2-32) 


(p  =  0 

n_  =  0 

n  •  V/?+  =  0 

n  '  f  e  =  '  f + 


On  the  anode,  all  the  positively  charged  ions  are  repelled.  The  negatively  charged  ions  and 
electron  number  densities  are  unaltered  when  reached  the  anode  surface.  The  electric  potential  on 
the  anode  is  the  balance  of  the  electromotive  force,  EMF  and  the  loss  along  the  external  circuit. 

<p  =  EMF  -  IR 

n  ■  Vn_  =  0 

-.-o'  (2'33) 

n  ■  Vne  =  0 


where  the  symbol  /  is  the  total  electrical  current  and  R  is  the  resistance  of  the  circuit  including 
the  discharging  column. 

The  boundary  conditions  on  the  far li eld  boundaries  are  uniformly  imposed  by  the  vanishing 
outward  nonnal  gradient  approximation  to  describe  the  vanishing  electrical  field  potential; 


n  ■V<p  =  0 
n  ■  Vn  =  0 
n  ■  V/7+  =  0 
n  ■  Vne  =  0 


(2-34) 


2 . 8  Characteristics  of  DCD 


The  typical  DCD  computational  simulations  are  conducting  at  the  ambient  pressure  levels  from 
nearly  vacuum  up  to  50  Torr  by  an  electromotive  force  (EMF)  up  to  9  kV.  The  electric  current 
density  is  usually  limited  around  10  mA/cm2.  The  classic  DCD  discharge  structure  has  been 
examined  by  the  parallel  electrodes  arrangement  and  is  implemented  by  the  side-by-side 
configuration  for  flow  control.  In  general,  the  dimension  of  discharge  domain  is  proportional  to 
the  magnitude  of  the  applied  EMF  and  inversely  proportional  to  the  ambient  pressure,  to  be 
precise,  the  ambient  gas  density.  The  important  current  voltage  relationship  has  been  studied  by 
the  classic  work  by  Engel  and  Steenbeck  in  one-dimension  [55]  and  most  recently  by  Rafatov  et 
al.  by  an  elaborate  drift-diffusion  model  [56].  At  the  low  current  density  and  due  to  the  greater 
electron  mobility;  the  ion  number  density  usually  exceeds  the  electron  counterpart  within  the 
entire  discharge  domain. 

Figure  3  depicts  the  contours  of  electron  and  positively  charged  number  density  of  the  DCD 
within  a  parallel  electrode  configuration.  The  diffusive  discharge  is  generated  between  a  gap 
distance  of  2.0  cm  by  an  EMF  of  2.4  kV  at  an  ambient  pressure  of  5  Torr.  On  the  cathode  all  the 
electrons  are  repelled,  but  the  secondary  emission  generates  a  rapid  avalanche  in  the  positive 
column.  The  positively  charged  ions  are  clustered  in  the  cathode  layer  to  reach  the  maximum 
number  density  of  3.0xl010/cm\  Although  the  electron  temperature  is  not  including  in  the 
simulation,  but  the  highest  value  from  the  classic  work  is  known  to  exist  in  the  cathode  layer  and 
has  a  value  around  3  eV  [7,55]. 
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Figure  3  Electron  and  Positively  Charged  Ion  Number  Density  Contours  of  Parallel  Electrodes 

p=5  Torr,  EMF=2.0  kV,  gap-2.0  cm 

In  Figure  4  by  increasing  the  ambient  pressure  from  3  to  10  Torr,  the  electric  field  intensity  is 
increasing  proportionally  [55].  The  distinct  cathode  fall  is  captured  by  the  two-dimensional 
computational  simulations.  In  fact,  the  intensities  drop  by  nearly  two  orders  of  magnitude  from 
the  cathode  layer  to  the  positive  column.  At  the  pressure  of  10  Torr,  the  maximum  normal 
electric  field  intensity  on  the  cathode  excesses  the  value  of  8. Ox  103  and  decreases  to  a  value  of 
3.5x10  V/cm  in  the  positive  column  and  rises  moderately  on  the  anode.  The  computational 
result  at  p=5  Torr  agrees  well  to  the  early  effort  by  a  different  numerical  algorithm  as  a 
validation  [4,5,33].  An  interesting  observation  is  also  made  in  that  by  lower  the  value  of  the 
secondary  emission  coefficient,  y  ;  the  thickness  of  the  cathode  layer  increases  accordingly.  This 
effect  also  produces  a  drop  in  electric  field  intensity  and  the  normal  current  density  which  means 
that  by  lowering  the  secondary  emission  coefficient  actually  reduces  the  ionizing  efficiency  in 
the  cathode  layer. 
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Figure  4  Surface  Normal  Electric  Field  Intensity 
EMF=2.0kV,  gap =2. Ocm 

In  the  DCD  modeling,  the  electrode  heating  is  rigorously  determined  by  the  Fourier’s  law  for 
conductive  heat  transfer,  and  the  dissipative  Joule  heating  can  be  evaluated  by  the  current  density 
and  the  electrical  field  intensity  induced  by  the  discharge.  The  temperature  of  the  neutral 
particles  has  been  simulated  by  Petrusiv  et  al.  [57]  by  assuming  a  constant  electron  temperature 
of  1  eV  (1.1604xl04  K).  It  is  anticipated  that  the  maximum  value  of  the  temperature  is  located 
deeply  in  the  cathode  layer.  Realistically,  the  nonequilibrium  electron  temperature  is  not  an 
invariant  in  the  discharge  domain  and  has  been  described  by  an  empirical  equation  of  reduced 
field  intensity  |  /P;  TJT  =  29.96  ln(|f?|  /  p)  +  24.64 .  The  typical  DCD  experiment  for  hypersonic 

flow  control  is  generated  around  an  EMF  of  2  kV  and  at  a  pressure  level  of  5  Torr  which  has 
exceeded  the  reduced  electrical  field  intensity  of  the  original  empirical  relationship.  For  the 
purpose  of  illustration,  the  electron  temperature  contour  of  a  DCD  is  depicted  in  Figure  5.  Again, 
the  maximum  electron  temperature  is  observed  to  be  embedded  within  the  cathode  layer  and 
decay  rapidly  in  the  positive  column  toward  the  anode  and  far  field.  The  range  of  electron 
temperature  variation  within  a  DCD  is  indeed  modest,  1  <  Te  <  2.8  eV. 
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X,  cm 


Figure  5  Electron  Temperature  Profile  with  a  DCD  Discharge 
EMF=2. 0  eV,  P=5.0  Torr 

The  DCD  structure  is  presented  in  Figure  6  by  a  side-by-side  electrodes  configuration  as  the 
most  flow  control  applications.  The  DCD  is  generated  at  an  EMF  of  2.4  kV,  a  pressure  level  of  5 
Torr,  and  with  the  gap  distance  between  electrodes  of  2.0  cm.  The  cathode  is  placed  on  the  left  of 
the  anode  in  the  composite  presentation.  On  the  left-hand  side  of  the  inset,  the  discharge  current 
vector  traces  are  appended  to  the  electron  number  density  contour.  On  the  right-hand  side  of  the 
inset,  the  electric  field  intensity  vector  traces  are  superimposed  over  the  positively  charged  ion 
number  density.  The  maximum  electron  density  over  the  anode  is  3.02x  10 10  and  the  positively 
charged  ion  is  located  over  the  cathode  with  a  maximum  value  of  5.76x  1010/cnf .  The  basic  DCD 
structure  is  still  identifiable  such  as  the  positive  column,  cathode  and  anode  layers  to  that  of  the 
parallel  electrodes  arrangement,  but  the  discharge  domain  becomes  pronounced  biased  by  the 
different  gap  distance  between  electrodes.  However,  the  average  cathode-layer  thickness  is  still 
maintained  at  0.125  cm,  which  is  comparable  to  the  result  of  the  parallel  electrodes  configuration. 
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Figure  6  Electron  and  Ion  Number  Density 

Contour  and  charge  current  density  (left),  ion  number  density  and  electric  field  intensity  (right); 

EMF=439V,  1=5.20  mA,  p=5  Torr 

It  is  also  interesting  to  note  that  the  strong  electric  field  vector  is  perpendicular  downward  to  the 
cathode  and  the  field  intensity  is  nearly  parallel  to  the  dielectric  surface  between  electrodes 
toward  the  cathode.  All  the  field  intensity  traces  are  highly  concentrated  at  the  edges  of  the 
electrodes  to  indicate  a  high  gradient  electric  field  region  locally. 

The  electron  number  density  obtained  from  Langmuir  probe  measurements  within  a  DCD  field  is 
shown  in  Figure  7.  A  60  mA  discharge  current  is  maintained  by  an  electrical  potential  at  1.2  kV 
over  the  side-by-side  electrode  arrangement  [33].  Surveys  are  conducted  on  the  plate  centerline 
and  the  cathode  is  placed  to  the  left  of  the  anode.  The  electron  number  density  is  presented  as  a 
function  of  x  at  four  different  heights  above  the  electrode  surface.  The  electron  number  density  is 
highest  near  the  plate  surface,  and  decreases  away  from  the  plate.  Peak  number  densities  of  about 
10  cm"  are  obtained  over  the  cathode,  and  the  measured  value  is  comparable  with  the 
numerical  results  under  the  similar  condition.  The  lowest  electron  concentrations  are  still 
measurable  at  0.5  cm  above  the  plate  surface  downstream  of  the  cathode,  but  the  number  density 
at  this  height  drops  off  rapidly  upstream  and  downstream  from  this  location 
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Figure  7  Measurement  of  Electron  Number  Density 
(over  a  side-by-side  DCD  configuration  by  a  Langmuir  Probe,  (p=I .2  kV,  1=60  mA,  p=1.0  Ton) 

The  behavior  of  the  electric  field  intensity  at  the  edges  and  over  the  electrodes  is  highlighted  by 
Figure  8.  The  two-dimensional  simulation  is  conducted  for  the  side-by-side  configuration  with 
the  electrodes  dimension  of  0.5  cm  in  length  and  the  gap  distance  of  1.0  cm.  The  outer  edges  of 
the  electrodes  are  placed  at  the  1.5  cm  from  the  outer  boundaries  of  the  computational  domain.  In 
order  to  resolve  the  high  intensity  field,  a  nonunifonn  250x120  mesh  system  is  used  with  a 
minimum  horizontal  grid  spacing  of  5><10"  cm  and  a  minimum  perpendicular  spacing  of 
1.671x0°  cm  clustering  to  the  electrodes.  The  normal  electric  field  intensities  over  the  cathode  at 
different  ambient  pressure  have  the  averaged  value  around  5kV/cm,  but  rise  sharply  to  four  times 
the  value  at  the  edges  of  electrode.  Similar  behaviors  are  also  exhibit  over  anode  but  with  a  much 
less  intensity.  This  phenomenon  is  well-known  in  electric  circuit  design  and  especially  in 
discharge  physics  [7,8]. 
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Figure  8  Normal  Electric  Field  Intensity 
EMF=2.0  kV,  gap=0.5  cm,  electrode  length=1.0  cm 

The  most  important  flow  control  mechanism  of  DCD  by  Joule  heating  distribution  is  depicted  in 
Figure  9.  The  volumetric  heating  occurs  mainly  within  the  cathode  layer  at  the  inner  edge  of  the 
cathode  toward  the  anode.  A  hot  spot  also  appeared  on  the  inner  edge  of  the  anode  but  with  a 
much  less  intensity.  From  the  numerical  simulation,  the  Joule  heating  release  to  the  discharging 
domain  is  1.3  J/s-cm2.  This  value  is  relatively  small  in  comparison  with  the  total  amount  of 
power  for  the  surface  DC  plasma  generation  of  21.52  J/s-ctn2.  The  Joule  heating  contributes 
about  50%  to  the  total  heat  flux  from  the  DCD  in  comparison  with  the  electrode  heating  at  the 
simulated  condition,  and  this  ratio  has  been  verified  by  experimental  measurements  [4,5,32].  The 
advantage  of  the  Joule  heating  is  beyond  the  electrode  surface  but  extends  into  the  flow  field.  In 
all,  the  net  effect  of  Joule  heating  is  created  a  local  thermal  expansion  to  the  flowfield  and  alters 
the  boundary-layer  displacement  thickness  and  thus  the  local  slope  of  a  control  surface.  The 
outward  displaced  surface  slope  always  triggers  a  family  of  compression  waves  by  the  external 
supersonic  stream  to  generate  the  flow  control  mechanism  [4,5,19]. 


26 

Approved  for  public  release;  distribution  unlimited 


1.5 


Y 

0.0 


QEJ 

8.30E-07 

4.68E-07 

2.57E-07 

1.42E-07 

7.79E-06 

4.28E-06 

3.36E46 

1.30E-06 

7.14E-05 

3.93E-05 

3.16E-05 

UK-05 

6.54E-04 


1.5  2  2.5  3 


X 


Figure  9  Joule  Heating  Contour 

(over  side-by-side  direct  current  discharge,  E=439.4  V,  1=5.20  mA,  p=5.0  Torr,  Gap=1.0cm.  Joule 
hearting  1.3  Watts,  DCD  power  supply  21.52  Watts) 

The  global  effect  of  the  externally  applied  magnetics  flux  density  to  DCD  and  the  comparison 
with  experimental  observation  is  given  by  Figure  10.  The  computational  simulations  is 
comparing  with  the  experimental  observation  in  at  Mach  number  of  5. 15  at  the  identical 
condition  of  a  pressure  of  1.0  Torr,  a  temperature  of  43  K,  and  by  an  EMF  of  1.2  KV.  The 
transverse  magnetic  flux  density  of  ±0.1  Tesla  is  applied  with  opposite  polarities  across  the  side- 
by-side  electrodes  by  a  gap  distance  of  2.45  cm.  At  the  relatively  low  externally  applied 
magnetic  field,  the  Lorentz  force  exerts  a  significant  influence  to  the  discharge  structure. 
According  to  the  right-hand  rule  and  in  the  specific  experimental  arrangement,  the  positive 
polarity  in  the  applied  transverse  magnetic  field  expels  the  plasma  away  from  the  electrode 
surface  by  the  Lorentz  force.  By  reverse  the  polarity  of  the  magnetic  field,  the  plasma  is 
suppressed  toward  the  electrode.  At  the  relatively  low  magnetic  flux  density  of  0. 1  Tesla,  the 
externally  applied  magnetic  field  exerts  a  profound  change  to  the  DCD  structure. 


Figure  10  DCD  in  Transverse  Magnetic  Field,  Bz=0,  Bz=-0.1,  and  Bz=+0.1 
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The  detailed  electrical  current  vectors  of  DCD  in  the  presence  of  an  externally  applied  transverse 
magnetic  field  are  given  in  Figure  11.  The  computational  simulation  mimics  the  experiment, 
except  by  a  smaller  computational  dimension  (electrodes  of  0.5  cm  in  length  and  1.0  cm  apart),  a 
weaker  magnetic  flux  density  of  0.05  Tesla,  but  at  a  stronger  EMF  of  2.0  kV.  The  corresponding 
Hall  parameter  of  the  electron  and  ion  are  p  =  /t  B  =  ±0.440  and  p.  =  ju.b .  =  ±0.0014 
respectively.  The  numerical  simulations  are  in  agreement  with  experimental  observation  in  that 
the  polarity  of  the  applied  magnetic  field  may  increase  plasma  generation  to  sustain  a  fixed 
discharge  electric  current.  In  the  numerical  simulation  at  B=-0.05  Tesla,  the  Lorentz  acceleration 
pushes  the  charged  particles,  thus  the  discharge  current  away  from  the  electrodes,  and 
concentrates  the  current  toward  the  inner  edges  of  the  electrode  pairs.  The  opposite  effect  is 
displayed  at  B=+0.05  Tesla:  The  discharge  electrical  current  is  suppressed  over  the  dielectric 
surface  between  electrodes.  The  discharge  current  is  nearly  parallel  and  confines  to  the  surface 
strait.  The  effect  for  a  greater  flow  control  effectiveness  has  also  been  recorded  by  experimental 
observations  [4,5,23-25,33]. 


Bz  =  -0.05  Tesla 

Bz  =  0.05  Tesla 

Bz  =  +0.05  Tesla 

Figure  1 1  Electrical  Current  Vectors  in  the  Presence  of  an  Externally  Applied  Magnetic  Field 
1.0  cm  separation  distance  between  electrodes,  EMF=2.0  kV,  P=5  Torr. 

2.9  Boundary  Conditions  of  Dielectric  Barrier  Discharge 

The  dielectric  coating  is  the  key  component  for  the  proper  DBD  operation,  because  it  limits  the 
amount  of  charge  transported  by  a  single  microdischarge  and  distributes  the  microdischarges 
over  the  entire  electrode  [9-11].  Therefore  one  of  the  electrodes  is  always  encapsulated  by 
dielectric  material  and  is  often  grounded  [3].  Despite  numerous  computational  simulations,  some 
computational  simulations  have  not  imposed  the  correct  media  interface  boundary  condition  for 
describing  this  key  physical  phenomenon  of  DBD  [44-46]. 

According  to  the  theory  of  electromagnetics,  the  electric  field  across  the  dielectric  and  plasma 
interface  must  satisfy  the  following  conditions  [19,50]; 

nx(Ed-Ep)  =  0,  (2-35) 

n-(Dd-Dp)  =  peil.  (2-36) 

In  the  above  equations,  the  subscripts  p  and  d  designating  the  variables  either  reside  in  plasma 
medium  or  on  the  dielectrics.  The  local  surface  charge  density  is  evaluated  by  integrating  the  net 
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charged  number  density  over  an  infinitesimal  normal  distance  over  the  control  surface,  and  is 
defined  as  pes  =  ej  (n  -  ne  -n_)dc>  which  has  a  dimension  of  coulombs/cm2. 

Equation  (2-35)  simply  states  that  the  tangential  electrical  field  strength  is  continuously  across 
the  media  interface.  It  can  be  shown  on  the  Cartesian  frame,  the  rate  of  change  for  electric 
potential  in  z,  dcp/dz,  must  be  identical  across  the  interface  along  the  x  coordinate.  Similarly,  the 
rate  of  change  of  electrical  potential  in  x,  d(p/d c,  must  be  equal  along  the  z  coordinate  across  the 
media  interface.  This  requirement  is  automatic  satisfied  by  the  two-dimensional  formulation. 

The  discontinuity  of  the  normal  component  of  the  electric  displacement,  D,  across  different 
media  must  be  balanced  by  the  net  surface  charge  on  the  interface  by  emission,  desorption,  and 
accumulation.  In  fact,  the  condition  defined  by  equation  (2-36)  is  independent  of  all  chemical- 
physics  processes  on  the  medium  interface.  And  this  equation  is  further  developed  by 
introducing  the  electric  potential  for  the  partially  ionized  plasma,  E  =  -V <p  to  become; 


d(PP  _  d(pd 


=  Pe 


(2-37) 


p  dn  dn 

In  the  above  equation,  £p  and  ed  denote  the  electric  permittivity  of  the  plasma  and  dielectrics 

respectively.  For  the  weakly  ionized  gas,  the  relative  value  of  £p  is  around  unity  and  for  has  a 
magnitude  up  to  7.0  from  polystyrene  to  rubber  [46], 


On  the  dielectric  barrier,  the  surface  diffusion  and  desorption  can  be  modeled  by  the  collisions  of 
excited  molecules  and  atoms  by  their  collision  frequency  and  the  binding  energy  [15].  The 
charge  accumulation  on  surface  is  considered  as  the  result  from  the  instantaneous  recombination 
of  the  charge  particles  after  satisfied  the  imposed  boundary  condition  for  the  species 
concentration.  For  the  two-dimensional  formulation,  equation  (2-37)  reduces  to  the  following 
fonn; 


d(PP  £d  d<Pd  ,  e  r, . w_ 

- = - 1 - (77 ,  —  n  —  n  )cl(7. 

dn  sp  dn  sp 


(2-38) 


Equation  (2-38)  is  the  critical  boundary  condition  to  be  imposed  on  the  dielectrics.  All  other 
physics  meaningful  boundary  conditions  on  the  interface  have  been  demonstrated  to  be  robust 
and  accurately  described  the  behavior  of  weakly  ionized  gas  generated  by  electron  collision  [43- 
46]. 


In  DBD  operation,  the  distinguish  discharge  phenomena  during  an  AC  cycle  is  frequently 
referred  to  as  the  forward  stroke  or  negative  going  and  the  back  stroke  or  positive  going  [10,1 1]. 
In  fact,  the  discharge  is  exclusively  governed  by  the  polarity  of  the  electrodes  [9,44-46].  When 
the  electrical  potential  is  positive,  (p(t)  >  0  on  the  exposed  electrode,  it  perfonns  as  the  anode  on 
which  the  electrical  potential  is  dictated  by  the  externally  applied  AC  current  source  with  a 
frequency  of  vj,  and  all  positively  charged  ions  are  repulsed.  The  imposed  boundary  condition  on 
the  anode  or  the  exposed  electrode  are; 

(p{t)  =  EMF  sin(s77),  (2-39) 

n+  =  0,  (2-40) 


29 

Approved  for  public  release;  distribution  unlimited 


n  ■  V«e  =  0  and  n  ■  V«_  =  0. 


(2-41) 


The  EMF  is  the  external  voltage  applied  to  the  DBD  circuit.  It  shall  be  immediately  recognized 
that  the  electric  field  potential  across  the  exposed  and  the  dielectric  barrier  after  the  breakdown  is 
determined  by  the  external  circuit  equation;  EMF  =  cp  +  R^  Jdv  [1 1,14,15],  and  the  electric 
current  density  of  the  discharge  is  simply,  J  =  e(f  _  f  -  f  _). 

At  this  same  instant,  the  dielectric  barrier  surface  acts  as  the  cathode  on  which  the  secondary 
electron  emission  is  induced  by  an  excessive  accumulation  of  positively  charged  ions  over  the 
dielectric  barrier  to  reduce  the  effective  difference  in  electric  potential.  At  this  AC  phase,  the 
imposed  boundary  conditions  on  the  cathode  for  the  electrical  potential  and  charged  particle 
concentration  on  the  dielectric  shall  be: 


n  ■  (V  (p)p  =  n  •  (fXY  <p)d  +  j-  J  (n+  -  ne  -  n_  )d  cr 


n  -  V«  =  0  and  n  =  0, 


(2-42) 

(2-43) 


f,=-rA-  (2-44) 

During  this  AC  cycle,  Equation  (2-42)  ensures  diminishing  electrical  field  intensity  across  the 
electrodes  through  the  surface  charge  accumulation.  It  is  also  the  fundamental  mechanism  of  the 
self-limiting  characteristic  in  preventing  DBD  transition  to  spark.  Again,  equation  (2-44) 
enforces  the  electron  secondary  emission  condition  for  the  Townsend’s  discharge  [7,8] 

When  AC  field  switches  polarity  to  a  negative  electric  potential,  cp(t)  <  0,  on  the  exposed 
electrode:  the  roles  of  the  exposed  electrode  and  the  dielectric  barrier  surface  also  reverses.  The 
exposed  electrode  now  functions  as  the  cathode;  the  secondary  emission  by  the  positive-ion 
accumulation  now  occurs  over  the  metallic  surface.  The  coefficient  of  the  emission  on  the  metal 
surface,  ym  is  different  from  that  of  the  dielectrics. 

During  this  AC  phase,  the  boundary  conditions  on  the  cathode  or  the  exposed  electrode  are; 

cp(t)  =  —EMF  sin  (cot)  (2-45) 

n  ■  Vn_  =  0  and  n+  =  0,  (2-46) 

r,  = -r.fi  •  (2-47) 


On  the  dielectric  surface  or  the  grounded  anode  at  the  same  time,  the  boundary  conditions  are; 
(Pit)  =  0,  (2-48) 

n+  =  o,  (2-49) 


n  ■  =  0  and  n  ■  V«_  =  0.  (2-50) 

To  complete  the  boundary  conditions  specification  for  the  DBD  simulation,  the  vanishing 
gradient  condition  is  imposed  unifonnly  on  the  farfield  boundary  of  the  computational  domain 
for  all  four  dependent  variables.  For  the  multiple  length-  and  time-scales  phenomenon;  the 
ionization  and  recombination  of  charge  particles  is  considered  to  be  instantaneous  in  comparison 
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with  the  time  scale  of  the  alternative  electric  field.  Therefore  no  separate  boundary  condition 
treatment  is  imposed  during  the  switching  AC  cycle. 

It  may  be  noticed  the  imposed  boundary  conditions  for  DBD  on  the  exposed  electrode  and  the 
dielectric  surface  are  fully  reinforced  the  DBD  self-limiting  mechanisms.  The  boundary 
conditions  specification  must  be  synchronized  with  the  alternating  polarities  of  the  electrode. 
The  differences  appear  only  in  the  secondary  emission  coefficients  on  the  cathode  and  the 
effective  electrical  potential  across  the  discharging  domain. 

2.10  Characteristics  of  DBD 


The  partially  ionized  gas  existing  in  the  discharging  domain  is  generated  through  a  succession  of 
random  streamer  fonnations  in  space  and  time  for  the  duration  of  few  nanoseconds  [9-11].  In 
application,  the  DBD  generally  operates  at  atmospheric  pressure  with  a  gap  distance  on  the  order 
of  one  to  a  few  millimeters  between  electrodes.  An  alternative  voltage  is  required  to  support  the 
random  transient  streamer  formation  in  the  electrode  gap  and  quenching  by  the  localized  charge 
build-up  on  the  dielectric  layers.  The  partially  ionized  gas  between  the  gap  of  electrodes  is 
generated  through  the  succession  of  a  large  volume  of  glow  discharge  and  random  streamers  in 
space  and  time  [3,9-1 1].  Gherardi  et  al.  [58]  recently  have  shown  that  when  the  DBD  occurs  in 
many  thin  filaments  will  lead  to  multiple  microdischarges  and  by  a  single  discharge  channel  over 
the  entire  electrode  surface  to  become  a  glow  discharge.  Through  the  emission  spectroscopy  and 
electrical  measurements,  the  transition  from  the  filamentary  structure  is  controlled  by  the  density 
of  metastable  nitrogen  molecules,  but  the  seeding  electrons  are  always  created  by  the  Penning 
effect. 

Specifically,  when  the  exposed  electrode  is  positively  biased  in  an  AC  cycle,  the  discharge  is 
characterized  by  a  streamer  like  structure.  In  the  negatively  biased  phase,  the  discharge  appears 
as  a  more  diffusive  structure.  In  essence,  the  DBD  discharges  are  operated  in  the  micro¬ 
discharge  mode  whose  nanoseconds  life  span  is  governed  by  the  ionization-recombination 
process.  On  the  other  hand,  the  propagation  of  charged  particles  is  synchronized  to  the  frequency 
of  the  AC  field.  Meanwhile  the  induced  electrostatic  force  by  the  free-space  charge  separation  in 
the  AC  cycle  becomes  a  repetitive  dynamic  event.  In  Figure  12,  this  dynamic  event  is  described 
in  a  composite  presentation  based  on  experimental  observations  in  a  typical  voltage-current 
curve  by  Moreau  [3]  and  high-speed  photography  by  Enloe  et  al.  [10,1 1]. 


Moreau,  grounded  encapsulated  electrode 


Figure  12  Voltage -current  Curve  and  Discharge  Photos  Exposed  Electrode  Anode  (U)  and  Cathode  (L) 
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The  important  comparison  of  the  periodic  voltage-current  characteristic  of  computational  results 
with  experimental  data  over  1,000  microseconds  (five  AC  cycles  at  5  kHz)  is  depicted  in  Figures 
13.  The  comparison  with  data  at  the  externally  applied  electrode  potential  of  4.0  kV  reveals  a 
good  overall  affinity  between  data  and  computational  results  at  the  instant  of  breakdowns.  Prior 
to  the  breakdown,  the  electrical  current  consists  only  of  the  displacement  component.  After  the 
plasma  ignition,  the  electrical  current  has  an  additional  conductive  component  of  the  DBD.  The 
magnitude  of  this  component  is  relatively  minuscule  and  indicates  a  current  density  less  than 
0.008  mA/cm.  This  calculated  value  agrees  well  with  data  collected  over  a  wide  group  of 
experiments  to  indicate  the  ionization  model  by  the  chemical-physics  kinetics  is  credible  [3,59]. 
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Figure  13  Compare  Voltage-current  Curves  (5kHz)  with  Experimental  Observation 

In  Figure  14,  the  fundamental  self-limiting  characteristic  of  the  DBD  is  convincingly  captured  by 
the  computational  simulation  at  the  externally  applied  electrical  potential  of  3.0  kV  and  5.0  kV 
across  the  overlapped  electrodes.  At  the  quiescent  atmospheric  condition,  the  discharge 
breakdown  takes  place  at  a  voltage  of  2.77  and  2.65  kV,  for  5.0  kV  and  3.0kV,  respectively, 
before  the  externally  applied  electrical  field  reaches  its  peak  values  for  both  polarities.  Then  a 
lower  and  constant  electric  potential  of  2.5  kV  is  maintained  by  the  conductive  current  and  by 
the  surface  charge  accumulation  on  the  electrodes.  The  numerical  result  duplicates  the  pattern  of 
experimental  surface  potential  measurements  [10,1 1].  In  the  positive  polarity  phase,  the 
discharge  occurs  when  the  positive-going  potential  exceeded  the  breakdown  voltage  and 
continues  until  the  negative-going  external  field  falls  beneath  the  breakdown  threshold.  The 
identical  behavior  is  also  observed  for  the  negative  polarity  phase.  During  the  initial  breakdown 
process,  a  sudden  drop  of  the  electric  potential  with  respect  to  time  induces  a  surge  of  the 
displacement  current  in  the  discharge.  In  experimental  measurements,  this  surge  indicates  the 
existence  of  multiple  micro  discharges  and  also  knows  to  induce  a  train  of  pulsations  in  the 
electrical  circuit  [3,9-1 1,59].  Similarly  the  sudden  drop  of  the  electric  potential  has  also 
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generated  non-physical  oscillations  by  the  numerical  simulation  known  as  the  Gibbs 
phenomenon  [37-41,43-45], 


Figure  14  Effect  of  Discharge  duration  with  Different  EMF 

The  domain  of  the  positively  charged  ion  at  the  peak  positive  voltage  is  compared  at  3kV,  5kV, 
and  8  kV  in  Figure  15.  The  discharge  domain  retains  the  similar  shape  but  its  dimension  is 
noticeably  affected  by  the  driving  voltage.  The  charged  ion  density  drops  exponentially  from  its 
peak  value  near  the  corner  of  the  exposed  electrode  and  the  dielectric  wall  to  the  ambient.  The 
maximum  positively  charged  ion  number  density  increases  from  4.3x10  ,  1.9x10  ,  to  4.0X  10  ~ 
per  cm3  with  the  increasing  intensity  of  the  applied  voltage.  The  number  of  electron  density  is 
also  increasing  proportional  with  the  rising  intensity,  but  at  a  two-order  of  magnitude  lower 
value  than  that  of  the  ion.  The  combination  of  a  higher  charge  number  density  and  a  longer 
discharge  duration  leads  to  a  great  periodic  electrostatic  force  to  agree  with  experimental 
observation  [3].  The  maximum  instantaneous  magnitude  of  the  force  generated  at  8  kV  is  more 
than  double  that  of  5kV,  to  be  as  high  as  2.7X105  dyne/cm3. 
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Figure  15  Discharge  Domains  of  DBD  by  Different  Applied  EMFs 


Figure  16  summarizes  the  effects  of  the  gap  distances  between  the  exposed  electrode  and  the 
encapsulated  electrode.  The  figure  shows  the  charge  density  and  electric  field  potential  for  the 
gap  distance  of  -1,  0  and  1  mm.  It  is  found  that  as  the  distance  increases  and  the  domain  of  DBD 
decrease,  and  the  effect  is  mainly  caused  by  the  changes  in  the  electric  potential  field. 
Overlapping  the  two  electrodes  can  slightly  enhance  the  DBD  effects  but  it  does  not  significantly 
increase  after  overlapping  the  two  electrodes  by  1mm.  These  observations  are  identical  for  both 
the  positive  and  negative  polarities  in  the  AC  cycle  and  agree  with  experimental  data  [59].  In 
essence,  the  computational  results  based  on  the  electrodynamics  honor  Paschen’s  law  [7].  One 
also  observes  that  the  discharge  is  persistently  for  either  the  positive-  or  negative-going  potential 
along  the  voltage-current  trace  after  the  breakdown  whether  in  the  positively  or  the  negatively 
biased  EMF  polarities. 
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Figure  16  Effect  of  Gap  Distances  to  DBD  Discharge 

From  the  experimental  observations,  the  relative  dielectric  permittivity  for  dielectric  materials 
has  exerted  measurable  influence  to  the  characteristic  of  discharge  [3,7,8,10].  The  dependence  of 
discharge  properties  to  the  dielectric  material  is  studied  by  the  model  equations.  Calculations 
were  perfonned  at  the  external  applied  voltage  of  4  kV  for  three  different  values  of  the  relative 
electric  permittivity,  s;  2.7  (polystyrene),  4.7  (glass),  and  6.7  (rubber).  Increasing  the  dielectric 
permittivity  leads  to  a  lower  threshold  of  breakdown  voltage  and  greater  amount  of  conductive 
current  of  the  discharge.  The  breakdown  voltage  drops  from  the  value  of  3. 1  kV  for  e  =  2.7  to 
2.7kV  for  e  =  4.7  and  then  to  2.6  kV  for  e  =  4.7.  As  can  be  seen  from  the  contour  plots  for  the 
positive  ion,  one  does  not  find  significantly  change  beyond  e  =  4.7.  The  lower  breakdown 
voltage  at  the  higher  values  of  the  electric  permittivity  results  in  an  earlier  ignition  and  a  delayed 
extinguishment;  thus  extends  the  discharge  duration.  The  higher  electric  permittivity  also 
produces  a  more  compact  discharge  domain  as  shown  in  Figure  17. 
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Figure  17  Effect  of  Permittivity  on  DBD  Domain 

The  final  and  the  most  important  verification  of  the  chemical-physics  ionizing  model  is 
presented  by  a  direct  comparison  with  a  measured  conductive  current  during  the  discharge.  The 
ionization  model  is  built  on  the  Townsend’s  discharge  mechanism  together  with  the  bulk  and 
ion-ion  recombination,  as  well  as,  the  electron  attachment  and  detachment.  If  the  ionization 
model  by  the  chemical-physics  kinetics  would  not  able  to  describe  a  reasonable  discharge 
composition,  the  calculated  conductive  current  of  a  DBD  will  be  clearly  erroneous.  The 
experimental  data  is  collected  at  an  AC  cycle  of  5  kHz  and  an  electrical  potential  of  4.0  kV  [59]. 
In  Figure  18,  the  calculated  conductive  current  is  designated  by  black  line  which  lies  within  the 
scattering  traces  of  the  experimental  data  to  reflect  the  correct  magnitude  and  duration  of  the 
discharge.  Equally  important,  the  effects  of  different  electrical  permittivity  and  secondary 
emission  coefficients  on  metallic  electrode  and  dielectric  surface  are  also  accurately  captured. 

4000 


2000 

> 

o 

0 

CD 

O 

> 

-2000 


-4000 

0  0.2  0.4  0.6  0.8  1 

Cycle  (ill) 

Figure  18  Comparison  of  Electrical  Conductivity  Current  with  Experimental  Measurement,  cp=4  kV, 

w=10 kHz 
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2.1 1  Periodic  Electrostatic  Force  of  DBD 

The  most  important  phenomenon  of  the  DBD  for  flow  control  is  the  periodic  electrostatic  force 
during  discharge.  A  range  of  formulations  have  been  put  forth  to  calculate  the  periodic  force  by 
Corke  et  al.  [1],  Boeuf  et  al.  [38-40],  Likanskii  et  al.  [42],  Shang  et  al.  [44-46],  Signh  [47],  and 
Jayaraman  et  al.  [60],  However,  the  dominant  mechanism  is  the  electrostatic  force  for  the 
aerodynamics-electromagnetics  interaction  [6]; 

F  =  peE  (2-51) 

From  the  charge  number  density  of  the  separated  free-space  charges,  the  period  force  can  be 
calculated  as: 


F  =  e(ne  +  n_  -  n+)V  (p 


(2-52) 


In  the  above  equation,  the  electrostatic  force  vanishes  either  in  the  charge-neutral  farfield  or 
when  the  gradient  of  an  applied  electrical  field  ceased.  The  maximum  instantaneous  magnitude 
of  the  force  based  in  the  normal  operational  range  of  DBD  is  calculated  on  the  order  of  104 
dyne/cm  at  an  EMF  of  4  kV  which  have  a  wide  range  of  variations  among  all  computational 
simulations  [38-40,  44-46].  However,  the  space-time  averaged  values  of  the  electrostatic  force 
yield  a  value  range  from  0  to  340  dyne/cm  in  the  positive  push  region  to  agree  well  with  the 
estimates  by  Beouf  et  al.  [38,39],  in  which  have  a  recorded  range  of  values  from  10  to  100 
dyne/cnr .  In  their  work,  even  an  instantaneous  peak  value  up  to  10  dyne/cnr  has  been 
estimated.  The  calculated  force  has  also  been  documented  in  the  physical  unit  per  electrode 
length  for  two-dimensional  computational  simulation,  and  has  a  value  up  to  24.4  dyne/cm  [40]. 
This  numerical  result  is  in  general  agreement  with  experimental  data  from  1.0  to  25  dyne/cm  by 
Hoskinson  et  al.  [61,62],  However  in  their  experiences,  the  exposed  electrode  is  constructed  by 
cylinders  of  varying  diameters. 

The  orientation  of  the  force  is  completely  detennined  by  the  gradient  vector  of  the  electric  field 
potential  and  the  sign  of  the  net  balance  of  charged  particles,  Equations  (2-5 1)  and  (2-52).  In  the 
two  dimensional  fonnulation,  the  force  components  in  the  x  and  y  coordinates  becomes; 


Fx=e(ne  +  n_-n+) 
Fy=e(ne  +  n_-n+) 


d(p 

dx 

d(p 

¥ 


(2-53) 


From  the  above  equation,  the  orientation  and  magnitude  of  the  resultant  periodic  electrostatic 
force  can  be  unambiguously  detennined.  The  uncertainty  can  only  be  challenged  by  the  physical 
fidelity  of  the  ionization  process  and  the  transport  properties  of  the  charged  gas  mixture.  The 
shortfall  in  basic  knowledge  of  the  chemical  kinetics  and  quantum  chemical-physical  kinetics  is 
the  stumbling  block  for  all  computational  simulations  [37-47]. 

The  Equation  (2-53)  shows  that  the  force  vector  must  follow  the  direction  of  the  electric  field 
vector;  it  is  useful  to  study  how  the  electric  field  intensity  is  modified  by  the  separated  free-space 
charges,  i.e.  the  net  balance  of  the  charged  number  density.  If  one  neglects  the  effects  of  the 
charged  density  and  normalized  the  electric  potential  field  by  the  instantaneous  potential  of  the 
exposed  electrode,  the  solution  is  dictated  by  the  externally  applied  AC  field,  as  shown  by  the 
subset  (a)  in  Figure  19.  The  contour  and  vector  plot  in  the  subset  (a-1)  are  the  dimensionless 
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electric  potential^,  and  electric  field  E  respectively,  and  contours  of  the  magnitude  of  the 
electric  field,  are  displayed  in  the  subset  (a-2).  The  subsets  (b)  and  (c)  in  Figure  19  show  the 
same  dimensionless  quantities  when  the  EMF  of  a  5  kHz  AC  cycles  are  in  its  maximum  and 
minimum  peaks  of  8.0  kV.  As  can  be  seen  from  the  figure,  the  separated  charge  number  density 
in  this  case  is  too  small  to  affect  the  solution.  The  instantaneous  force  field  is  acting  along  a 
fixed  externally  applied  electric  field  potential  and  alters  direction  when  the  polarity  on  the 
exposed  electrode  changes  sign. 


Figure  19.  Comparison  of  Dimensionless  Electric  Potential 
Electric  field  and  its  magnitude  for  (a)  no  changed  density,  (b)  maximum  penitential  and  (c)  minimum 

potential. 

Attention  is  now  focused  on  the  separated  charged  density  term  in  the  RHS  of  Equation  (1-53). 

In  here,  the  discussion  is  limited  to  the  nL.-n+  model.  A  more  complete  discussion  including  n 
will  be  elaborated  later  but  the  conclusion  is  identical.  Figure  20  depicts  the  comparison  of 
number  density  of  the  positively  charged  ion  and  electrons  in  the  peak  positive  and  negative  AC 
cycle.  In  the  positive  cycle,  the  exposed  electrode  is  the  anode,  shown  in  subsets  (a-1)  to  (c-1); 
the  positively  charged  ions  are  repelled  by  the  exposed  electrode  and  are  attached  to  the 
dielectric  surface.  The  positively  charged  ions  propagate  along  the  electric  field  lines  and  are 
accumulated  on  the  dielectric  surface.  The  distributions  of  positively  charged  ion  density  (red 
color)  and  electron  density  (black  color)  are  also  inserted  in  the  subset  (b-1).  The  magnitude  of 
positively  charged  ion  density  is  always  two  orders  of  magnitude  greater  than  that  of  the  electron 
number  density  at  the  peak  AC  cycle.  This  indicates  that  the  separated  charge  number  density, 
and  hence  the  force,  is  dominated  by  the  positively  charged  ion  number  density.  In  subset  (c-1), 
the  force  vector  follows  the  electric  field  intensity  line  and  its  magnitude  indicated  by  the 
contour  that  resembles  the  shape  of  the  positively  charged  ion  density  contour.  Similar 
observation  can  be  made  in  the  peak  negatively  biased  AC  cycle,  as  shown  in  subset  (c-1)  to  (c- 
2). 
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Figure  20.  Comparison  of  (a)  Positive  Ion  Density,  (b)  Electron  Density,  and  (c)  Force  Field 
The  insert  in  (b-1)  shows  the  maximum  value  for  the  positive  ion  density  (red)  and  electron  density 

(black). 


The  global  structure  of  the  electric  field  intensity  is  easily  understood  because  when  the  exposed 
electrode  is  anode,  the  electric  field  intensity  decreases  toward  the  dielectrics;  thus,  the  electric 

field  intensity  is  positive  over  the  dielectrics  (E  —  — V tp,E> 0).  The  electric  field  vector  over  the 
dielectric  must  synchronize  with  the  changing  AC  polarity,  and  becomes  negative  ( E  <  0  )  when 
polarity  reversed.  Therefore,  the  relative  number  densities  of  electrons,  positively  and  negatively 
ions  ultimately  determine  the  orientation  of  the  force.  This  observation  has  been  verified  over 
the  entire  range  of  the  externally  applied  electrical  potential  examined,  at  different  gap  distances 
between  electrodes,  and  variations  of  electrical  permittivity  of  dielectrics  [43-46].  Therefore,  the 
orientation  of  the  instantaneous  and  the  time-averaged  periodic  electrostatic  forces  is  no  more 
complex  than  just  the  push-pull  or  pull-push  combination.  It  solely  depends  on  the  net  balance  of 
the  positively  and  negatively  charged  number  density  and  the  local  electric  field  intensity. 

Another  interesting  issue  is  the  effect  of  the  negatively  charged  ions  to  the  periodic  electrostatic 
force.  To  address  this  question,  one  needs  to  examine  the  generation  and  depletion  processes  of 
the  negatively  charged  ions.  The  theory  of  nonequilibrium  chemical  kinetics  indicates  the  key 
mechanisms  are  the  electron  attachment,  detachment,  and  ion-ion  recombination.  From  most 
detailed  analyses  [13-15],  the  number  density  of  the  negatively  charged  ion  of  DBD  is  much 
lower  in  comparison  with  the  positively  charged  counterpart  -  generally  appears  as  a  fraction  of 
a  few  percent.  Again  all  computational  simulation  within  the  EMF  ranging  from  3.0  to  8.0  kV, 
the  number  densities  of  the  negatively  charged  ions  are  nearly  two  orders  of  magnitude  lower 
than  the  positively  charged  counterpart;  the  typical  values  are  around  n+  a  ~  o(io12 cm-3)  and 

n  max  ~  o(\Q'"Cin  3)  •  These  values  agree  very  well  with  the  most  elaborated  calculations  from  the 
chemical  kinetic  models  by  Bogdanov  et  al.  [13]. 
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In  Figure  21,  this  quantified  difference  in  charged  particle  number  densities  difference  has  been 
further  confinned  by  the  electrical  current-voltage  curve  during  an  AC  cycle  of  5  kHz  and  driven 
by  an  EMF  of  8  kV.  The  computational  results  by  the  two-  and  three-component  ionization 
models  are  undiscemible.  To  demonstrate  the  maximum  possible  effects  of  negative  ion,  the 

solution  with  kd  =10  14  is  depicted.  The  electrical  field  intensity  shows  an  identical  breakdown 
voltage  of  2.9  kV.  During  the  discharge  phase  for  both  the  three-  and  two-species  ionization 
models  indicate  the  same  constant  voltage,  2.4kV,  the  more  sensitive  conductive  current  also 
shows  identical  behavior  for  the  two  ionization  models,  which  means  the  negatively  charged  ions 
contribute  little  to  the  discharge  current. 


Cycle  (t/T) 

Figure  21.  Electric  Current  Voltage  Patterns  of  the  Two-charged-species  (ne,  n+)  and  Three-charged- 
species  (ne,  n_,  n+)  Ionization  Model,  EMF=8.0  kV,  f=5.0  kHz 

Figure  22  displays  the  comparison  between  the  two-  and  three-species  ionization  models  at  the 
instance  the  electric  potential  over  the  exposed  electrode  is  positive  and  reaches  its  maximum 
value.  The  data  are  depicted  in  a  composite  presentation  for  the  electric  field  intensity  vector  E, 
the  net  balance  of  different  number  density  of  charged  particles  (ne  -  n_)  and  (n+  —  n_  —  ne),  as 
well  as,  the  instantaneous  electrostatic  force  vector  at  an  applied  EMF  of  8.0  kV  across  the 
overlapping  electrodes.  First  of  all,  the  discharge  patterns  are  unaltered;  all  the  electric  field 
intensity  vectors  point  along  the  electric  field  line,  shown  in  (a-1)  and  (b-1).  The  instantaneous 
and  maximum  value  of  (n+  -  ne)  over  the  dielectric  barrier  attains  a  value  of  3.96x  1012  cm'3.  On 
the  other  hand,  the  (n+  —  ne  —  n.)  indicated  a  value  of  3.97x  10  “  cm'  which  reflects  a  relatively 
insignificant  negatively  charge  ion  number  density.  The  resultant  force  according  to  equation 
(8-3)  is  unifonnly  directed  from  the  exposed  electrode  to  the  dielectric  barrier  and  has  an 
instantaneous  maximum  magnitude  of  2. 133x10  dyne/cm  .  The  numerical  result  is  nearly 
identical  to  the  previous  result  of  2.200x  10  dyne/cm  by  a  two  species  charged  modeling  [45]. 
This  result  is  greater  than  the  maximum  force  of  4.10xl04  dyne/cm3  generated  by  the  EMF  at  4 
kV,  and  of  1.500x10'’  dyne/cnr1  by  the  EMF  of  6.0  kV  [31,45].  However,  it  needs  to  be 


40 

Approved  for  public  release;  distribution  unlimited 


reminded  that  the  maximum  electrodynamics  force  only  exists  within  the  thin  plasma  sheath  and 
for  a  very  short  period  during  the  discharge. 


Figure  22.  Comparison  of  Two-  and  Three-species  Ionization  Models  of  DBD  at  the  Peak  Positive  Cycle 
(a)  two-species  model  and  (b)  three-species  model.  Subset  1  is  the  electrical  field  lines  and  contour  is  the 
electron;  subset  2  is  the  negative  ion  density  contour.  The  insert  in  subset  2  is  the  comparison  of  the 
maximum  values  of  electron  density  (red)  and  negative  ion  density  (black);  subset  3  is  the  positive  ion 
contour.  The  insert  in  subset  2  is  the  comparison  of  the  maximum  values  of positive  ions  (blue)  and  the 
sum  of  electron  and  negative  ion  densities  (black);  subset  4  contain  the  Electrodynamic  force  vectors  and 

charge  density  contour.  V=8.0  kV,  f=5.0kHz 

Figure  23  displays  the  electric  potential  and  charged  particle  number  density  distributions  for  the 
instantaneous  DBD  field  when  the  exposed  electrode  is  acted  as  the  cathode.  All  ions  are 
repulsed  from  the  dielectric  surface  and  clustered  around  the  lower  corner  of  the  exposed 
electrode.  At  the  same  instant,  a  large  number  of  faster  moving  electrons  are  propagating  further 
over  the  dielectric  barrier  away  from  the  comer  region  of  the  electrodes.  The  sign  of  electric  field 
intensity  is  now  negative;  pointed  from  the  dielectric  barrier  toward  the  exposed  electrode. 
According  to  the  present  computation,  the  plasma  sheath  is  extremely  thin  (less  than  0.002  mm) 
and  its  thickness  decreases  with  a  lower  value  of  the  electron  secondary  emission  coefficient. 

The  maximum  value  of  charge  separation,  (n+  -  ne)  near  the  intersection  of  electrodes  is 
4.4018x10  cm'  ,  the  corresponding  net  number  density  of  charged  separation  of  the  three- 
species  model  (n+-  n_  -  ne)  is  4.4059x10  “.  As  a  consequence,  the  forces  generated  by  the  two 
models  are  very  also  very  close,  1.0012xl03  dyne/cm3  and  9.9879xl04  dyne/cnr  for  the  three 
and  two  species  models,  respectively.  At  this  instance  of  the  AC  cycle,  the  peak  instantaneous 
electrostatic  force  is  greater  than  when  the  electrical  potential  on  the  exposed  electrode  is 
positively  biased  but  the  maximum  force  occurs  on  the  surface  of  the  exposed  electrode. 
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Figure  23  Comparison  of  Two-  and  Three-species  Ionization  Models  of  DBD  at  the  Peak  Negative  Cycle 
(a)  two-species  model  and  (b)  three-species  model.  Subset  1  is  the  electrical  field  lines  and  contour  is  the 
electron;  subset  2  is  the  negative  ion  density  contour.  The  insert  in  subset  2  is  the  comparison  of  the 
maximum  values  of  electron  density  (red)  and  negative  ion  density  (black);  subset  3  is  the  positive  ion 
contour.  The  insert  in  subset  2  is  the  comparison  of  the  maximum  values  of positive  ions  (blue)  and  the 
sum  of  electron  and  negative  ion  densities  (black);  subset  4  contain  the  Electrodynamic  force  vectors  and 

charge  density  contour.  V=8.0  kV,  f=5.0  kHz 

When  the  exposed  metallic  electrode  acts  as  the  cathode,  the  periodic  electrostatic  force  is 
oriented  from  the  dielectrics  toward  the  exposed  electrode.  The  magnitude  of  this  force  is 
generated  by  a  sufficiently  high  electrical  intensity  and  charged  particles  number  density  within 
the  intersection  of  electrodes.  This  magnitude  of  the  electrodynamics  force  is  not  negligible,  but 
it  is  constrained  by  the  solid  surface  of  electrodes  to  become  ineffective  for  momentum  transfer 
with  the  neutral  species.  Unfortunately,  all  the  key  elements  are  exclusively  dependent  on  the 
chemical  and  chemical-physical  kinetics  processes  that  are  not  fully  understood  to  become  the 
key  inaccurate  contributor  to  DBD  computational  simulations  at  the  present. 

The  time-averaged  periodic  electrostatic  forces  by  the  two-  and  three-charged  species  ionization 
models  over  an  AC  cycle  and  the  entire  discharge  domain  is  depicted  in  Figure  24.  The  overall 
force  structure  is  indistinguishable  from  the  two  ionization  models  due  to  the  effect  of  the 
number  density  of  negative  ions  is  not  significant  as  compared  to  that  of  the  positive  ion.  In  short, 
by  considering  the  presence  of  the  negatively  charged  ions  with  a  chemical-physics  kinetics 
model,  the  averaged  force  that  produces  the  well-known  electric  wind  is  less  than  without  it,  but 
is  also  insignificant.  However,  it  must  be  emphasized  that  the  precise  quantification  still 
uncertainty  must  be  accepted  as  the  state-of-the-art  for  DBD  simulations. 
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Figure  24  Time -averaged  Electromagnetic  Force  of  DBD  over  a  Completely  AC  Cycle,  EMF=8.0  kV, 

f=5.0kHz 

There  are  two  distinct  regions  in  the  direction  of  the  time-averaged  forces  separated  by  a  null- 
value  line  cutting  diagonally  across  the  comer  from  the  exposed  electrode  and  the  dielectric  wall. 
The  region  below  the  line  of  demarcation  is  the  positive  force  creating  a  mean  pushed  motion 
along  the  downstream  direction.  The  region  above  the  cutting  line  has  the  negative  force  and  its 
maximum  force  is  exerting  toward  the  vertical  surface  of  the  exposed  electrode. 

2.12  Concluding  Remarks 

The  drift-diffusion  model  for  plasma  flow  control  is  a  viable  approximation  for  the  transport 
properties  of  low-temperature  partially  ionized  gas  by  electron  collisions.  All  ionizing 
phenomena  are  occurred  across  the  molecular/atomic  scales  and  beyond  the  scope  of  gas  kinetics 
theory;  thus  must  be  modeled  either  by  chemical-physics  kinetics  or  nonequilibrium  chemical 
kinetics.  Despite  this  glaring  deficiency  for  plasma  modeling  and  simulation,  the  fundamental 
formulation  still  can  be  used  to  delineate  the  ambiguities  from  both  the  experimental 
observations  and  computational  simulations. 

In  DCD  application,  the  dominant  process  of  the  Joule  and  electrode  heating  and  the  Joule 
heating  is  concentrated  over  the  cathode  and  mostly  from  the  inner  edge  of  the  side-by-side 
cathode-anode  configuration.  The  heat  released  by  the  Joule  heating  is  generally  lower  than  six 
percent  of  the  power  input  for  DCD  generation. 

The  multiple  microdischarges  within  the  DBD  domain  are  unsupportable  by  current 
computational  capability.  As  the  consequence,  the  different  DBD  mechanisms  during  an  AC 
cycle  are  unable  to  be  duplicated  precisely  by  the  drift-diffusion  and  chemical-physics  models. 
However,  the  periodic  electrostatic  force  in  different  phases  of  the  AC  cycle  still  can  be 
described  by  the  polarities  of  the  exposed  electrode  with  the  grounded  and  encapsulated 
electrode. 
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The  confusion  issue  on  the  orientation  of  the  periodic  force  can  be  clarified  by  the  net  balance  of 
the  number  density  of  electrons,  the  negatively  and  positively  charged  ions.  The  instantaneous 
periodic  electrostatic  force  is  no  more  complex  than  by  a  simple  push-pull  or  push-push 
description.  The  maximum  time-averaged  force  has  a  lower  magnitude  than  the  instantaneous 
counterpart.  According  to  the  chemical-physics  ionization  models,  the  negatively  charged  ion 
contributes  very  little  to  the  global  DBD  perfonnance. 

A  sustained  basic  research  is  urgently  needed  for  future  advance  of  nonequilibrium  chemical 
reaction  and  quantum  chemical-physical  kinetics  at  full  spectrum  of  thermal  conditions. 
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3  Modeling  Fuel  Cells 


3.1  Introduction 

In  order  to  investigate  SOFCs  mathematically,  efforts  have  been  put  into  development  of  models 
including  mass  transportation  and  electrochemical  reactions.  Zhu  et  al.  [63]  presents  a  new 
computational  modeling  framework  for  SOFC  simulation  that  takes  the  whole  system  including 
flow  channels  and  planar  membrane-electrode  assemblies  into  consideration.  His  work  employed 
multistep  reaction  mechanisms  in  terms  of  detailed  elementary  heterogeneous  chemical  kinetics. 
Detailed  charge  transfer  reactions  are  analyzed  by  separating  the  mechanism  into  several 
elementary  steps.  Greene  et  al.  [64]  focuses  on  minimized  the  concentration  overpotential  by 
applying  functionally  graded  electrodes  and  observes  the  physical  phenomenon  of  mass  transfer 
throughout  the  electrodes  for  multi-gas  inputs.  However,  the  activation  overpotential  in  both  of 
their  models  is  directly  calculated  by  the  Butler- Volmer  equation  and  do  not  take  the  micro 
structural  characteristics  into  consideration.  Theoretically,  these  micro  structural  factors  are 
critical  to  the  size  of  active  reaction  surface  sites  and  hence  affect  the  rate  of  electrochemical 
reaction.  M.  Ni  et  al.  [65]  developed  the  model  from  a  micro-scale  level  and  the  model  was  able 
to  capture  the  coupled  electrochemical  reactions  and  mass  transfer  involved  in  SOFC  operation. 
But  all  the  micro  structural  parameters  are  treated  separately.  Although  physically  all  those 
parameters  are  observed  to  correlate  with  each  other.  Once  one  parameter  changes,  the  rest  of  the 
parameters  should  alter  correspondingly.  This  work  expands  upon  previously  developed  theories 
and  models  [66,67].  The  model  takes  into  account  electronic,  ionic,  and  gas  transport  together 
with  electrochemical  reaction  effects.  It  can  predict  the  distribution  of  overpotentials,  current 
densities,  and  gas  concentrations  along  the  electrode.  Also,  the  model  takes  into  account  all  the 
microstructural  factors  that  are  critical  to  cell  performance.  In  addition,  the  model  applied  the 
binary  random  packed  sphere  model  to  mimic  the  microstructural  make-up  the  electrodes  [68,69]. 
However,  to  the  best  knowledge,  the  literatures  lack  exploring  the  correlation  of  microstructural 
parameters,  which  is  crucial  to  resemble  a  real-world  cell  performance. 

The  primary  focus  of  our  study  is  to  investigate  how  the  microstructural  parameters  are  related  to 
each  other  and  how  they  affect  cell  perfonnance.  Several  studies  associated  with  tortuosity  and 
porosity  relations  were  developed  and  organized  by  Matyka  et  al.  [71].  For  a  spherical  particle 
mixture,  Currie  [70]  proposed  that  porosity  is  inversely  proportional  to  tortuosity.  This 
assumption  was  confirmed  by  experimental  data.  Ricardo  Dias  et  al.  [71]  extended  this  study  and 
investigated  the  adjustable  parameter,  n,  which  is  crucial  to  determining  the  porosity-tortuosity 
correlation  (see  equation  3 1).  German  [73]  summarized  that  porosity  is  dependent  on  the  relative 
composition  of  two  species  and  the  particle  size  ratio  in  a  binary  mixture  of  spherical  particles. 
The  larger  the  particle  size  ratio  is  (i.e.  the  larger  the  difference  between  the  sizes  of  the  two 
particles),  the  higher  the  packing  density  at  all  compositions  will  be.  Ricardo  Dias  et  al.  [74] 
explored  the  dependence  of  packing  porosity  on  particle  size  ratio,  focusing  on  exploring  the 
significance  of  particle  size  ratio. 

3.2  Numerical  Model 

3.2.1  Anode 

3.2. 1 . 1  Oveipotential  due  to  electrochemical  reactions  and  ohmic  resistance 
The  overall  charge  balance  relationship  can  be  written  as: 
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(3-1) 
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Where  Je  a  and  Jt  a  are  the  current  density  (A/m“)  due  to  transport  of  electronic  and  ionic 

2  3 

conductors  in  the  anode;  Sv  is  active  surface  area  per  unit  volume  (m  /m  )  of  the  porous 
electrode;  J  n  a  is  the  transfer  current  density  per  unit  area  of  reaction  surface  (A/m  ).  This 
equation  accounts  for  the  electrochemical  reaction  rate  for  the  fuel  cell  along  the  anode. 

The  active  surface  area  per  unit  volume  indicates  the  available  reactions  sites  that  can  be  used  for 
electrochemical  reactions  and  was  developed  by  Costamagna  et  al.  [69]  from  binary  random 
packing  and  percolation  theories.  It  is  represented  below: 


7  7 

S„  —  n  sin2  Or^nnn ,  — PP 


(3-2) 


where  0  is  the  contact  angle  between  electronic  and  ionic  conducting  particles;  r  is  the  radius  of 
electronic  conducting  particles;  nt  is  the  total  number  of  particles  per  unit  volume;  «,•  and  ne  are 
the  number  fraction  of  ionic  and  electronic  conducting  particles,  respectively;  Z.  and  Zg  are  the 

coordination  number  of  ionic  and  electronic  conducting  particles,  respectively;  Pt  and  Pe  are 
the  probability  that  a  given  ionic  or  electronic  conducting  particle,  respectively,  belongs  to  a 
percolation  cluster.  Essentially,  Sv  can  be  expressed  as  a  function  of  porosity,  electronic  particle 
size,  number  fraction,  and  particle  size  ratio  [69]. 

The  transfer  current  density  is  normally  described  by  the  general  fonn  of  the  Butler-Volmer  (B- 
V)  equation. 
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/•  ■  2  *  * 
o  a  is  the  exchange  current  density  of  the  anode  (A/nr);  yH  is  the  molar  fraction  of  //, 

y'H^  is  the  molar  fraction  of  at  the  fuel  channel.  R  is  the  ideal  gas  constant  (8.3 14  J/mol  K). 

T  is  the  operating  temperature  (K).  (3  is  the  charge  transfer  coefficient  and  is  normally  chosen 
to  be  0.5  for  symmetry  [75]. 

Applying  Ohm’s  law  for  the  electronic  and  ionic  conductors,  we  get: 
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pfa  is  the  effective  resistivity  (Dm)  of  the  anode  electronic  conductors;  pfa  is  the  effective 

resistivity  of  anode  ionic  conductors;  Ve  and  Vt  are  the  electronic  and  ionic  potential  (V), 
respectively.  The  effective  resistivity  can  be  determined  as: 


peff  = 

r  e,a 


-■pf  = 

’  Hi. a 


(3-5) 


Where  is  the  volume  fraction  of  electronic  conductors;  T  is  tortuosity  of  the  anode;  S  is 
porosity  of  the  anode;  cre  a  is  the  electronic  conductivity  (S/m)  of  the  anode  electronic  conductors; 

and  <Jj  a  is  ionic  conductivity  (S/m)  of  the  anode  ionic  conductors. 

The  anode  overpotential  ,j  can  be  determined  by  the  difference  of  equilibrium  potential  and 
practical  potential. 

n.=ty? -v?)-(v,-v)  (3-6 ) 


C  and  Vr  are  the  equilibrium  electronic  and  ionic  potential  (V),  respectively,  and  are 
constant  throughout  the  cell. 

The  first  derivative  of  lja  can  be  written  as: 
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Combing  the  charge  balance  and  B-V  equation,  the  second  derivative  of  na  is  equal  to: 
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3.2.2  Overpotential  due  to  Mass  transport 

Diffusion  processes  within  a  porous  electrode  structure  can  be  distinguished  as  two  types.  First, 
there  is  nonnal  diffusion  in  which  one  gas  diffuses  through  another,  with  negligible  influence  of 
the  pore  walls  on  the  rate  of  diffusion.  This  applies  when  the  mean  free  path  of  the  molecules  is 
much  less  than  the  pore  diameter.  Second,  when  the  mean  free  path  of  the  molecules  is  greater 
than  the  pore  diameter,  this  is  referred  to  as  Knudsen  diffusion.  For  most  of  SOFCs,  the  Knudsen 
effects  cannot  be  neglected  [76].  Therefore,  for  a  binary  gas  system  going  through  the  pore 
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structure,  the  overall  effective  diffusion  coefficient  Df  can  be  written  by  combing  effective 
onnal  binary  diffusion  coefficient  Df  and  the  effective  Knudsen  diffusion  coefficient  D'f  .  [77] 
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(3-9) 


The  effective  diffusion  coefficient  depends  on  the  microstructure  of  the  porous  anode,  quantified 
through  the  porosity  and  tortuosity  values.  Thus,  the  effective  binary  diffusion  coefficient  can  be 
written  as: 
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For  Knudsen  diffusion,  the  coefficient  can  be  described  as  in  [63], 
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where  d0  is  the  pore  diameter  (m)  and  is  assumed  to  be  approximately  equal  to  the  hydraulic 
diameter.  [78] 
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A{)  is  the  specific  surface  area  based  on  the  solid  volume.  For  randomly  packed  binary  spheres, 
A0  is  expressed  as: 


=  6  ne  +  (l  -ne)a 
de  ne+(l-ne)a~3 


where  de  is  the  diameter  (m)  of  electrically  conducting  particles,  CL  is  the  particle  size  ratio  of 
ionic  to  electronic  conducting  particles. 

Similarly  as  with  the  effective  binary  diffusion  coefficient,  the  effective  Knudsen  diffusion 
coefficient  can  be  expressed  as: 
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The  general  fonn  of  the  Fick’s  law  takes  into  account  both  diffusion  and  convection  mass 
transfer  and  can  be  expressed  as  [79]: 
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(3-15) 
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Where  N.  is  molar  flux  (mol/m2  s)  of  species  i;  Df  is  the  effective  diffusion  coefficient  (m2/s) 

•5 

of  species  i  (from  equation  10);  Ct  is  the  concentration  (mol/nr)  of  specie  i;  v  is  the  convection 
velocity  (m/s). 

Under  constant  operating  temperatures,  the  ideal  gas  law  can  be  written  as: 


dp  _  dc  RJ, 
dx  dx 


(3-17) 


If  pressure  is  unifonn  throughout  the  electrode,  then  dp/dx  is  constant.  According  to  eqn.  (3-17), 
dc/dx  will  also  be  constant.  Then  we  will  have: 


dcH 


dx  dx 


(3-18) 


For  equimolar  counter-current  mass  transfer,  we  have  NHi  =  - Nn  o .  Combine  this  condition  with 
eqn.  (3-15)  and  eqn.  (3-16),  the  FF  flux  can  be  expressed  as: 


Nu 


(3-19) 


Based  on  flux-current  relations  as  well  as  ideal  gas  law,  we  have 
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Eqn.  (3-20)  turns  into: 
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3.2.3  Anode  governing  equations  and  boundary  conditions 

Combining  the  equations  and  solving  for  the  molar  fraction  of  fuel,  electronic  current  density, 
and  overall  overpotential;  we  get  a  total  of  three  governing  equations  for  the  anode. 
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The  computational  domain  is  shown  in  Figure  25.  The  boundary  conditions  for  the  governing 
equations  can  be  derived  as  follows:  at  the  fuel  gas  inlet,  which  is  also  the  location  of  the  current 
collector,  the  hydrogen  molar  fraction  is  equal  to  the  bulk  flow  value.  The  total  current  density 
only  comes  from  the  transport  of  electrons.  As  a  result,  the  boundary  conditions  can  be  expressed 
as: 
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Figure  25  Anode  Computational  Domain 

At  the  electrode-electrolyte  (EE)  interface,  the  transport  of  ions  is  the  only  factor  that  contributes 
to  the  overall  current  density.  Therefore,  ion  current  density  equals  to  the  overall  current  density. 
This  defines  the  boundary  condition  of  the  overpotential  at  the  EE  interface. 
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After  solving  the  three  coupled  governing  equations,  t]a  >  J e,a  and  >’//,  distributions  can  be 
obtained.  Now  the  overall  overpotential  of  the  anode  can  be  written  as  follows: 
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3.3  Cathode 

3.3.1  Overpotential  due  to  electrochemical  reactions  and  ohmic  resistance 

The  electrochemical  reaction  equations  in  the  cathode  are  similar  to  the  ones  in  the  anode  and 
can  be  derived  in  a  similar  fashion,  to  produce  the  following: 
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where  T)c  is  the  cathode  overpotential,  J 0  c  is  the  cathode  exchange  current  density,  pfc  is  the 

effective  resistivity  of  cathode  electrically  conducting  particles,  and  Je  c  is  the  electronic  current 
density. 

3.3.2  Overpotential  due  to  mass  transport 

Following  the  derivation  in  Berger  [80],  an  effective  Knudsen  diffusion  coefficient  of  O2  can  be 
defined  as: 
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The  effective  normal  diffusion  of  oxygen  taking  both  conduction  and  convection  transport  can  be 
defined  by 
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s  normal 


Total  concentration  of  O2  is  equal  to  the  summation  of  Knudsen  concentration  and  normal 
concentration. 
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On  the  cathode  side,  nitrogen  does  not  involve  any  electrochemical  reaction,  so  flux  of  nitrogen 
at  steady  state  is  zero.  Then  we  have  Ntotal  =  N ()_  .  The  O2  flux  turns  into 
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By  combining  flux-current  relations  and  the  ideal  gas  law,  Eqn.  (2-28)  becomes: 
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3.3.3  Cathode  governing  equations  and  boundary  conditions 

Similar  to  the  anode,  three  coupled  governing  equations  for  the  cathode  can  be  defined. 
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The  boundary  conditions  in  the  cathode  can  also  be  obtained  using  the  same  approach  used  with 
the  anode. 
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After  solving  the  governing  equations,  the  cathode  overall  overpotential  can  be  calculated  as: 

7C|  overall  =  (K  '  ~  K  ?)_  ^Fe\ inlet  ~K\ Ee)=  7C|  inlet  +  7C|  EE  (3-30) 
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3.4  Analysis  of  Microstructure  Parameters  Correlations  and  Percolation 
Threshold 

3.4. 1  Relationship  between  porosity  and  tortuosity 

Figure  26  plots  the  data  from  the  experiments  of  Currie  [70]  for  mixtures  of  spherical  particles. 
The  tortuosity  (T)  is  assumed  to  be  inversely  proportional  to  porosity  (£).  By  applying  eqn.  (2-31) 
as  the  relationship,  it  can  be  observed  from  Figure  26  that  most  of  the  data  falls  in  the  region 
where  the  n  value  is  between  0.4  and  0.5. 

r  =  4r  (3-31) 


e 


Figure  26  Dependence  of  Tortuosity  on  the  Packing  Porosity 

In  order  to  further  investigate  the  adjustable  parameter  n,  Ricardo  [71]  perfonned  an 
experimental  study  of  binary  mixtures  of  spherical  particles,  as  shown  in  Figure  27.  In  this 
figure,  cpL  represents  volume  fraction  of  large  particles.  It  can  be  observed  that  for  particle  size 
ratios  that  are  less  or  equal  to  3.33,  the  value  of  n  is  approximately  0.5  regardless  of  the  volume 
fraction.  For  most  SOFCs,  the  sizes  of  ionic  and  electronic  conductors  are  comparable,  so  it  is 
reasonable  to  choose  a  value  of  n  to  be  equal  to  0.5  in  the  numerical  model.  However,  if  we 
happen  to  encounter  a  large  particle  size  ratio,  the  plot  in  Figure  27  will  be  used  for  selecting  the 
appropriate  n  value. 
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Figure  27  Dependence  of  n  on  cpL  for  Different  Particle  Size  Ratios 
3.4.2  Relationship  between  porosity  and  particle  size 

For  common  SOFCs  the  particle  sizes  of  electronic  and  ionic  conductors  will  be  close,  and  hence 
the  ratio  will  be  close  to  1.  S.  Yerazunis  et  al.  [80]  performed  experiments  to  analyze  the  binary 
spherical  mixtures  with  comparable  particles  sizes  as  shown  in  Figure  28.  The  particles  were 
densely-packed  to  ensure  good  connection  of  particles,  as  in  most  electrode  configurations  of 
SOFCs.  Based  on  the  experimental  data,  porosity  can  be  expressed  as  a  function  of  particle  size 
ratio  and  volume  fraction  a  shown  in  eqn.  39. 
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Where  cp  is  volume  fraction  and  the  subscript  L  and  S  indicates  the  larger  particle  size  and 
smaller  particle  size,  respectively.  It  can  be  deduced  from  the  experimental  data  that  the  porosity 
value  approaches  0.36  as  the  particle  size  ratio  approaches  1.  Therefore,  we  assume  that  when 
ionic  particle  size  is  equal  to  electronic  particle  size,  the  porosity  is  set  to  be  0.36.  It  is  noted  that 
in  this  particular  circumstance  bimodal  mixtures  can  be  expressed  as  monomodal. 
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Figure  28  Packing  Density  vs.  Composition  for  Different  Particle  Size  Ratios 

3.5  Percolation  Threshold 

In  this  model,  the  electrode  was  assumed  to  be  represented  as  a  collection  of  randomly  packed 
spherical  particles  made  up  of  either  electronically  conductive  or  ionically  conductive  materials. 
By  applying  a  coordination  number  model  together  with  percolation  theory,  this  model 
guarantees  that  the  same  type  of  particles  (ionic  or  electronic  conductors)  contact  each  other  and 
form  a  network  or  particle  chain  through  the  electrode.  That  is  to  say,  for  any  given  combination 
of  compositions  and  particle  size  ratios  of  electronic  and  ionic  conductors,  the  coordination 
model  can  be  used  to  differentiate  if  it  is  above  or  below  the  percolation  threshold.  If  above  the 
percolation  threshold,  pathways  are  guaranteed  to  be  fonned  through  the  electrode  for  the 
conduction  of  ions  and  electrons.  If  below  the  percolation  threshold,  these  pathways  are  not 
guaranteed  to  span  the  entire  electrode  and  may  not  provide  continuous  pathways  for  conduction. 
In  these  cases,  the  model  assumptions  are  not  valid. 

3.5.1  Model  V alidation 

In  order  to  detennine  if  the  correlations  of  microstructure  parameters  have  the  potential  to 
accurately  predict  cell  performance,  model  validation  was  carried  out.  Three  experimental 
investigations  from  literature  [81]  [82]  [83]  focusing  on  different  perspectives  of  SOFC 
performance  were  selected  to  facilitate  the  comparison  between  predicted  results  and 
experimental  data. 

There  are  two  types  of  inputs  to  the  model:  operational  inputs  and  physical  inputs.  The 
operational  inputs  of  the  model  include:  operating  temperature  and  pressure;  fuel  gas 
composition;  and  current  density.  The  physical  inputs  include:  exchange  current  density; 
thicknesses  of  the  anode,  cathode  and  electrolyte;  volume  fraction  (or  mass  fraction)  of 
electronic  and  ionic  conductors;  particle  size  of  electronic  (or  ionic)  conductors;  particle  size 
ratio  of  ionic  to  electronic  conductors;  porosity;  and  tortuosity. 


55 

Approved  for  public  release;  distribution  unlimited 


The  first  experimental  data  selected  is  from  S.P.  Jiang  et  al.  [81].  The  objective  of  the  experiment 
was  to  investigate  the  effect  of  impregnation  of  different  volume  fractions  of  nano-sized  YSZ 
particles  into  nickel  anodes  on  the  electrode  behavior.  Table  2  lists  the  parameters  provided  by 
the  experimental  study.  In  order  to  solve  this  problem  using  the  developed  numerical  model,  Ni 
particle  size,  tortuosity,  and  anode  exchange  current  density  needed  to  be  determined. 


Table  1  e,  and  ^j/kB  values 


n2 

o2 

ch4 

h2o 

CO 

h2 

co2 

3.798 

3.467 

3.758 

2.641 

3.69 

2.827 

3.941 

^i/kB 

71.4 

106.7 

148.6 

809.1 

91.7 

59.7 

195.2 

Table  2  Value  of  input  parameters  for  model  validation  case  No.  1 


Parameter  (provided  by  paper) 

Value 

Operating  Temperature 

1073K 

Pressure 

1.0  atm 

Anode  Thickness 

30  pm 

Electrolyte  Thickness 

limn 

Anode  Gas  Composition 

97%  H2  (3%  H20) 

Volume  Fraction  of  NiO/YSZ 

Ni:  100%/0% 

Ni+2.7mg/cm2  YSZ:  83%/17% 
Ni+4.0mg/cm2  YSZ:  79%/21% 

Porosity 

30% 

YSZ  Particle  Size 

0. 1-0.3  pm 

2 

Let  us  take  the  Ni+4.0mg/cm“  YSZ  case,  for  example.  The  volume  fraction  of  Ni/YSZ  for 
Ni+4.0mg/cm2  YSZ  is  given  as  79%  Ni/21%  YSZ.  The  number  fraction  of  the  electronic 
conductor  can  be  expressed  as  a  function  of  particle  size  ratio  from  eqn.  (3-33).  The  probability 
that  one  cluster  that  belongs  to  a  percolation  cluster  can  be  obtained  from  eqn.  (3-34).  Then  the 
upper  and  lower  bounds  of  particle  size  ratio  and  electronic  number  fraction  are  calculated  and 
listed  in  Table  3.  The  values  within  the  two  bounds  will  maintain  percolation  and  hence  ensure 
good  conductivity. 
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n 


<Pe  = 


(1  -ne)a3  +ne 


P.  = 


(4-zY-5 

l  2  J 


1 0.4 


where 


(3-33) 


(3-34) 


z  6n‘- 

[n,+{\-n)a2 

Substituting  the  porosity  value  and  the  volume  fraction  of  the  electronic  conductor  into  the 

porosity-particle  size  ratio  correlation  gives  us  a  particle  size  ratio  of  ^=0.4.  By  choosing  the 
average  YSZ  particle  size,  which  is  0.2  pm,  Ni  particle  size  can  be  derived  as  0.5  pm. 

Table  3  Upper  and  lower  bounds  for  particle  size  ratio  and  electronic  volume  fraction  for  79%  Ni/  21% 


YSZ 


^=0.79 

Lower  bound 

Upper  bound 

a 

0.13 

0.53 

ne 

0.09 

0.36 

Tortuosity  can  be  obtained  from  the  porosity-tortuosity  correlation.  Since  the  particle  ratio  in  this 
case  is  less  than  3.33,  it  is  reasonable  to  set  n  equals  to  0.5. 


T 


1 

0.30'5 


1.826 


Exchange  current  density  was  calculated  from  the  provided  electrochemical  impedance 
spectroscopy  (EIS)  data  and  charge  transfer  resistance  is  determined  by  applying  the  simplified 
B-V  equation.  After  all  the  parameters  applied  in  the  model  are  calculated,  the  mathematical 
simulation  is  performed.  From  Figure  29,  it  can  be  seen  that  the  numerical  results  and 
experimental  data  agree  reasonably  well  with  each  other. 
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Figure  29  Comparison  between  Numerical  Results  and  Experimental  Data  for  Case  No.  1 

The  second  experiment  comes  from  Kim  et  al  [82].  The  objective  of  this  experiment  was  to 
investigate  the  performance  and  durability  of  Ni-coated  YSZ  anodes  for  intermediate 
temperature  solid  oxide  fuel  cells.  Applying  the  same  approach  as  the  previous  case,  ionic 
particle  size  in  anode,  electronic  particle  size  in  cathode,  and  tortuosity  is  calculated  by  using  the 
adopted  sub-model  correlations.  The  anode  and  cathode  exchange  current  density  are  calculated 
by  the  provided  EIS  plot  from  experimental  measurement.  After  that,  the  current-voltage  (I-V) 
curve  is  calculated  by  applying  the  mathematical  model.  Again,  a  good  consistency  between 
numerical  results  and  experimental  data  for  both  of  the  operating  temperatures  can  be  seen  from 
Figure  30. 


Figure  30  Comparison  between  numerical  results  and  experimental  data  for  case  No.  2 

The  last  experimental  data  is  from  S.  P.  Jiang  [83].  The  objectives  of  his  experiments  were  to 
examine  the  influence  of  sintering  temperature  on  cell  performance  by  using  an  anode  half-cell. 
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The  paper  provided  SEM  pictures  of  the  anode  microstructure  at  four  different  sintering 
temperatures.  In  Figure  31,  the  electronic  and  ionic  particle  sizes  are  approximated  from  the 
SEM  pictures.  The  green  and  red  circles  are  used  denote  Ni  and  8  mol%  Y203-ZrC>2  (YZ8Y) 
particles,  respectively.  The  range  of  particle  size  ratio  that  can  maintain  the  percolation  threshold 
is  calculated  to  be  between  0.5  and  2  for  this  volume  fraction  composition.  At  1,300  °C,  YSZ 
particle  size  is  approximately  0.5  pm,  while  the  Ni  particle  size  is  about  2pm  resulting  in  a 
particle  size  ratio  of  r  =  0.25.  Under  these  conditions,  the  assumptions  of  the  model  break  down 
because  the  percolation  threshold  is  not  satisfied.  Therefore  the  predicted  results  would  not  be 
valid.  However,  we  can  still  use  the  closest  limiting  particle  size  ratio  (i.e.  r  =  0.5)  to  calculate 
an  I-V  curve  to  compare  with  the  experimental  data.  For  this  case,  the  mathematical  predictions 
should  have  better  performance  than  experimental  data  because  percolation  threshold  is  assumed 
to  be  satisfied.  This  same  approach  is  used  for  the  1,350  °C  sintering  temperature  case,  as  well, 
because  of  the  extreme  particle  size  difference  of  two  types  of  conductors  (0.6  pm  for  Ni  and  2 
pm  for  YSZ  resulting  in  r  =  0.3).  For  1,400  °C,  the  Ni  and  YSZ  particle  size  are  estimated  as 
1.25  pm  and  1.5  pm,  respectively,  resulting  in  r  =  0.83.  For  1,500  °C,  both  Ni  and  YSZ  particle 
size  are  estimated  as  2  pm,  resulting  in  r  =  1.  The  particle  size  ratios  are  within  the  percolation 
threshold  for  the  1,400  and  1,500  °C  sintering  temperature  cases  and  the  numerical  simulation 
can  be  applied.  Once  the  particle  sizes  of  electronic  and  ionic  conductors  are  detennined,  the 
porosity  is  obtained  by  using  the  porosity-particle  size  ratio  correlation.  The  method  of 
calculating  tortuosity  in  this  case  is  a  little  different  from  the  previous  ones,  in  order  to  better 
match  the  provided  experimental  results.  Since  the  paper  provides  the  values  of  anode  ohmic 
resistance,  Equation  (5)  can  be  applied  to  estimate  the  tortuosity.  However,  if  the  porosity- 
tortuosity  correlation  is  still  applied,  the  estimated  tortuosity  would  be  an  order  of  magnitude  less 
the  values  which  correspond  to  the  resistance  measurements. 


59 

Approved  for  public  release;  distribution  unlimited 


Figure  31  SEM  of  Ni/  8  mol%  Y203-Zr02  (Ni/TZ8Y)  Cermet  Anodes  Sintered  at  (a)  1300°C,  (b)  1350°C, 

(c)  1400°C, (d)  1500°C 

Figure  32  shows  the  comparison  of  model  predictions  vs.  experimental  data.  It  can  be  seen  that 
for  1,400  and  1,500  °C  sintering  temperatures,  the  numerical  results  matched  reasonably  well 
with  the  experimental  data.  However,  for  the  remaining  two  sintering  temperatures,  there  is  a 
noticeable  difference  between  model  results  and  experimental  data,  as  expected.  For  these  cases, 
the  numerical  model  calculated  a  scenario  where  particle  size  ratio  was  artificially  selected  to 
satisfy  the  percolation  threshold.  However,  in  reality  the  particle  size  ratio  falls  outside  the 
percolation  threshold  range  indicating  poor  connectively  between  the  particles  of  the  ionic  and 
electronic  conductors.  Therefore,  observed  experimental  data,  where  the  terminal  voltage  drops 
rapidly  and  reaches  zero  at  a  very  small  current  density  value,  follows  with  expected  behavior 
based  on  our  model  assumptions. 
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Figure  32  Comparison  between  Numerical  Results  and  Experimental  Data  for  Case  No.  3 

3.6  Model  Sensitivity  Study 

The  model  validation  showed  a  reasonable  ability  to  predict  experimental  performance  results 
from  given  cell  microstructural  parameters.  Next,  a  sensitivity  study  is  perfonned  to  investigate 
the  effect  of  sub-model  correlations  on  model  predictions.  The  first  experimental  case  from  the 
model  validation  section  is  selected  as  a  baseline.  This  is  because  the  experiment  has  been  tested 
and  validated  indicating  a  good  consistency  with  the  results  from  the  numerical  model.  Also,  it 
involves  an  anode  half-cell  only,  so  there  is  no  need  to  worry  about  differentiating  effects 
between  the  anode  and  cathode.  Selected  variables  will  be  varied  to  test  the  sensitivity  of  model 
predictions  to  these  variations,  while  the  rest  of  the  parameters  will  be  obtained  from  the 
experiment  validation  case.  In  the  sensitivity  study,  two  different  volume  fractions  are  selected  to 
perfonn  the  tests  in  order  to  expand  the  sample  sizes. 
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able  4  Value  of  input  parameters  for  model  validation  case  No.  2 


Parameter  (provided  by 
paper) 

Value 

Operating  Temperature 

700°C  /  800°C 

Pressure 

1.0  atm 

Anode  Thickness 

1200pm 

Electrolyte  Thickness 

7ptn 

Cathode  Thickness 

30  pm 

Anode  Gas  Composition 

97%  H2  (3%  H20) 

Cathode  Gas  Composition 

air 

Volume  Fraction  of  NiO/YSZ 

40%/60% 

Mass  Fraction  of  LSM/YSZ 

50%/50% 

Porosity 

40% 

NiO  Particle  Size  (anode) 

20-30nm 

YSZ  Particle  Size  (anode) 

<300nm 

3.6.1  Tortuosity  vs.  Porosity 

There  are  two  places  containing  the  tortuosity  parameter  in  the  model.  They  are  related  to  ohmic 
resistance  in  electrode  and  mass  diffusion.  First  of  all,  since  electronic  conductivity  is  several 
orders  of  magnitude  higher  than  the  ionic  conductivity,  ohmic  overpotential  is  primarily 
dominated  by  the  resistance  due  to  ionic  conduction.  Based  on  Equation  (2-5),  it  can  be  found 
that  increasing  tortuosity  will  result  in  an  increase  in  effective  ionic  resistivity  and  hence  worsen 
the  perfonnance  of  the  SOFC.  Physically,  this  follows  from  the  fact  that  as  tortuosity  increases, 
the  path  length  for  conduction  from  one  end  of  the  electrode  to  the  other  end  increases.  Secondly, 
the  parameter  of  tortuosity  also  appears  in  the  effective  diffusivity  to  impact  the  gas  diffusion 
process.  From  the  eqn.  (2-10),  it  can  be  seen  that  the  increase  of  tortuosity  will  result  in  a 
decrease  in  the  effective  diffusion  coefficient  (or  a  corresponding  increase  in  mass  transfer 
resistance).  In  that  case,  a  larger  tortuosity  will  lead  to  a  lower  concentration  of  fuel  at  electrode¬ 
electrolyte  interface  and  hence  increase  the  concentration  overpotential. 
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Table  5  Value  of  input  parameters  : 

for  model  validation  case  No.  3 

Parameter  (provided  by 
paper) 

Value 

Operating  Temperature 

1273  K 

Pressure 

1.0  atm 

Anode  Thickness 

20-40pm 

Electrolyte  Thickness 

0.9±0.03mm 

Anode  Gas  Composition 

97%  H2  (3%  H20) 

Volume  Fraction  of  NiO/YSZ 

50%/50% 

NiO  Particle  Size 

Approximate  from  SEM 

YSZ  Particle  Size 

Approximate  from  SEM 

For  this  test  case,  the  porosity  is  given  as  30%.  The  idea  is  to  vary  the  unknown  parameter  n  in 
the  porosity-tortuosity  correlation  given  in  eqn.  (2-31)  in  order  to  find  out  how  tortuosity  impacts 
the  model.  The  n-value  is  chosen  to  vary  from  0.01  to  1  for  79%/21%  Ni/YSZ  and  from  0.01  to 
1.3  for  83%/17%  Ni/YSZ,  respectively.  The  overpotential  values  calculated  from  the 
corresponding  varying  n  values  are  compared  with  the  baseline  experimental  data,  which  is 
represented  by  the  model  with  n  =  0.5  (which  was  used  for  the  model  validation  case).  The 
purpose  is  to  explore  how  the  n-value  affects  the  predicted  anode  overpotential.  Figure  33  shows 
a  plot  of  tortuosity  (which  is  a  function  of  n)  versus  the  %-error  difference  when  compared  with 
the  baseline  experimental  data  (i.e.,  n  =  0.5).  The  dashed  region  indicates  the  range  of  ±5% 
difference  between  the  predicted  and  experimental  values  of  anode  overpotential.  In  order  to 
maintain  the  predicted  overpotential  between  the  5%  error  bars,  the  range  of  the  n-value  is 
approximately  from  0.275  to  0.685  for  79%/21%  Ni/YSZ  and  from  0.08  to  0.79  for  83%/17% 
Ni/YSZ,  respectively.  Correspondingly,  the  tortuosity  value  and  percentage  range  is  about  1.4  - 
2.25  (-23.3%  -  33.7%)  for  79%/21%  Ni/YSZ,  and  1.1  -  2.5  (-39.6%  -  36.9%)  for  83%/17% 
Ni/YSZ.  This  illustrates  that  a  relatively  large  variation  in  the  tortuosity  value  (~  +/-40%)  results 
in  a  small  impact  on  the  predicted  overpotential  (+/-  5%). 
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(b)  83%/17%Ni/YSZ 


Figure  33  Comparison  of  Predicted  Anode  Overpotential  at  Different  n-values 
3.6.2  Particle  size  ratio  vs.  porosity 


The  anode  governing  equations  indicates  that  both  mass  diffusion  and  electrochemical  reaction 
rate  depend  on  the  particle  size  ratio.  For  79%/21%  volume  fraction  Ni/YSZ,  the  range  of 
particle  size  ratio  is  calculated  to  be  between  0.13  and  0.53  to  maintain  percolation  threshold. 

For  83%/ 17%  volume  fraction  Ni/YSZ,  the  range  of  particle  size  ratio  is  calculated  to  be 
between  0. 1  and  0.4  to  maintain  percolation  threshold.  This  means  that  the  ionic  conducting 
particles  will  always  be  smaller  than  the  electronic  conducting  particles  for  this  case.  Based  on 
the  calculation  from  the  numerical  model,  it  can  be  found  that  as  the  particle  size  ratio  increases 
(getting  closer  to  1),  the  active  surface  area  decreases  and  results  in  an  increase  in  cell 
overpotential  as  shown  in  Figure  34.  This  is  because  as  particle  size  ratio  (a=ri/re)  gets  smaller 
(ionic  particles  get  smaller  or  electronic  particles  get  larger),  the  size  difference  of  two  types  of 
conductors  is  expected  to  be  more  significant.  Therefore  electrochemical  reactions  rate  will  be 
boosted  by  increased  contact  area  between  the  electronic  and  ionic  conducting  materials  at  the 
TPB  region  and  hence  improve  the  cell  performance.  For  mass  transfer,  as  particle  size  ratio 
increases  (getting  closer  to  1),  the  diffusion  coefficient  will  get  smaller  based  on  eqns.  (2-11),  (2- 
12),  and  (2-13),  leading  to  an  increase  in  the  concentration  overpotential.  However,  in  this  case, 
the  particle  size  ratio  impact  for  gas  diffusion  is  not  as  significant  due  to  an  extremely  thin  anode 
thickness. 


(a)  79%/2 1  %  Ni/YSZ  (b)  83%/l  7%  Ni/YSZ 


Figure  34  Effect  of  particle  size  ratio  on  active  surface  area 
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For  both  of  the  test  volume  fractions,  the  varying  factor  will  be  particle  size  ratio  only.  The 
particle  size  ratio  is  chosen  to  vary  from  0.35  to  0.44  for  79%/21%  Ni/YSZ  and  the  experimental 
data  is  represented  by  selecting  a=0.4,  which  was  used  for  the  model  validation  case.  For 
83%/ 1 7%  Ni/YSZ,  the  particle  size  ratio  is  varied  from  0.34  to  0.4  and  the  experimental  data  is 
represented  again  by  the  a  value  used  for  the  validation  case,  which  was  0.382.  As  shown  in 
Figure  35,  if  the  overpotential  difference  between  numerical  results  and  experimental  data  is 
expected  to  be  within  5%  error  (dashed  line  region),  the  deviation  of  particle  size  ratio  and  its 
percentage  needs  to  be  0.386  -  0.414  (-3.5%  -  3.5%)  for  79%/21%  Ni/YSZ,  and  0.376  -  0.385  (- 
1.3%  -  1.1%)  for  83%/17%  Ni/YSZ,  respectively.  For  both  of  the  two  volume  fractions,  the 
model  performance  predictions  are  strongly  sensitive  to  the  particle  size  ratio. 


Figure  35  Comparison  of  Predicted  Anode  Overpotential  at  Different  Particle  Size  Ratio  Values 


3.6.3  Discussion 

To  sum  up,  Table  6  lists  the  results  of  the  sensitivity  study  for  the  sub-model  correlations  at  two 
different  volume  fractions.  It  can  be  seen  that  the  model  is  extremely  sensitive  to  the  porosity- 
particle  size  ratio  correlation.  These  parameters  can  only  be  deviated  within  a  narrow  range, 
which  is  within  ±3.5%  in  order  to  maintain  a  ±5%  variation  in  the  predicted  overpotential  for 
both  of  the  volume  fraction  cases.  Therefore,  it  is  critical  that  this  value  to  accurately  detennined 
and  represented  in  the  model.  However,  the  model  does  not  seem  to  be  very  sensitive  to  the 
porosity- tortuosity  correlation.  For  a  small  variation  in  the  predicted  overpotential  (±5%),  the 
tortuosity  can  have  a  deviation  of  up  to  ±20%  to  ±30%.  Therefore,  the  porosity-tortuosity 
correlation  adapted  here  should  be  valid  to  approximate  tortuosity  for  model  predictions. 
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Table  6  Results  of  sensitivity  study 


Percent  change  corresponding  to  ±5%  change  in  predicted  overpotential  compared  with 
experimental  data 

Sample 

79%>Ni  /  21%  YSZ 

83%  Ni  / 17%  YSZ 

Tortuosity 

-23.3%  -  33.7% 

-39.6%  -36.9% 

Particle  size  ratio 

-3.5%  -  3.5% 

-1.3%  -  1.1% 

3.6.4  Conclusion 

In  this  study,  SOFC  performance  was  investigated  numerically  by  applying  two  sub-model 
correlations  (porosity-tortuosity,  porosity-particle  size  ratio)  aiming  to  tie  microstructural 
parameters  to  physically  measurable  parameters.  The  model  was  validated  against  available 
experimental  data  and  demonstrated  the  potential  to  predict  real  SOFC  perfonnance  numerically 
rather  than  spending  tremendous  time  on  performing  expensive  experiments.  The  sensitivity 
analysis  was  performed  on  the  adopted  sub-model  correlations.  It  was  determined  that  the 
porosity-particle  size  ratio  correlation  had  a  greater  effect  on  predicted  cell  performance  results 
than  the  porosity-tortuosity  correlation.  Since  the  model  has  the  capability  to  differentiate  which 
microstructural  parameter  plays  a  more  important  role  in  affecting  the  cell  performance,  this 
study  has  the  potential  to  help  improve  cell  efficiency  and  optimize  cell  performance  by 
adjusting  microstructural  parameters. 
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4  Fabrication  and  Testing  of  Advanced  Solid  Oxide  Fuel  Cells 


4. 1  Impact  of  In-line  Mixer  for  Printing  Anode  Interlayer: 

One  of  the  problems  with  printing  anode  interlayer  is  the  fact  that  the  two  phases,  NiO  and  YSZ 
do  not  mix  well  [84]  Layering  effects  were  seen  with  phase  segregation  of  NiO  and  YSZ.  In 
order  to  improve  the  mechanical  mixing,  a  mixer  insert  (a  Teflon  spiral)  was  introduced 
downstream.  Four  button  cells  were  fabricated,  two  of  which  served  as  control  cells.  The  anode 
interlayers  for  the  control  cells  (a)  and  (b)  were  printed  without  the  use  of  the  mixer  insert. 

Anode  interlayers  for  cells  (c)  and  (d)  were  printed  using  the  mixer  insert.  The  anode  interlayers 
were  printed  on  a  NiO-YSZ  support  substrate.  To  print,  NiO-YSZ  interlayers,  NiO  ink  and  YSZ 
inks  were  prepared.  YSZ  ink  was  prepared  by  mixing  the  powder  (Yttria  stabilized  zirconia)  with 
solvent,  plasticizers  and  dispersants.  NiO  ink  was  prepared  in  a  similar  manner  using  NiO 
powder.  Details  of  ink  preparation,  mixing  etc.  can  be  obtained  from  ref  (1).  The  same  procedure 
for  mixing  described  in  ref.  (1)  was  used.  Several  passes  of  the  anode  interlayer  was  printed  on 
the  NiO-YSZ  support.  For  the  electrolyte  printing,  YSZ  ink  was  prepared.  Several  passes  of  YSZ 
was  printed  on  the  anode  interlayer.  Traditional  slurry  pasted  LSM/YSZ  cathode  interlayer,  and 
cathode  current  collection  layer,  LSM  was  used  to  complete  the  cell  fabrication.  Out  of  the  four 
cells,  one  control  cell  failed.  Results  of  the  remaining  three  cells,  (a),  (c)  and  (d)  are  shown  in 
Figure  36  The  power  densities  of  the  cells  fall  in  the  range  from  350  mW/cm  to  430  mW/cm2  at 
850  °C.  Cells  (a)  and  (c)  show  very  similar  values  of  current  density  and  power  density.  The 
current  density  at  850  °C,  at  600  mV  overpotenial  is  about  800  rnA/cm2.  Cell  (d)  shows 
somewhat  improved  perfonnance.  The  current  density  at  850  °C,  at  600  mV  overpotenial  is 
about  1  A/cm“.  The  scatter  in  the  results  makes  it  difficult  to  make  a  conclusive  statement  about 
the  impact  of  the  in-line  mixer.  These  experiments  will  be  repeated.  In  order  to  evaluate  the 
microstructure  of  these  cells,  scanning  electron  microscopy  was  carried  out  after  electrochemical 
characterization.  Fig  4-2  shows  the  scanning  electron  image  (SEI)  of  the  cross-section  of  a 
typical  cell  with  anode  interlayer.  A  dense  electrolyte  can  be  seen,  about  12  pm  thick.  The 
remaining  layers,  anode  interlayer  Ni-YSZ,  and  cathode  layers,  LSM/YSZ  and  LSM  are  porous. 
The  cathode  interlayer  and  cathode  current  collection  layer  are  not  readily  distinguishable  from 
each  other.  Similarly  the  anode  interlayer  is  not  distinguishable  from  the  anode  substrate.  Fig  4-3 
shows  the  SEI  of  the  region  immediately  below  the  electrolyte.  In  the  case  of  cell  (a)  that 
contains  no  anode  interlayer,  the  region  seen  is  the  anode  substrate,  Ni-YSZ.  In  the  case  of  cells 
(c)  and  (d),  the  region  seen  is  the  anode  interlayer,  presumably  10-20  pm  thick.  The  region  is 
porous.  However,  the  Ni  and  YSZ  phases  are  not  distinct.  The  phase  contrast  requires  low 
voltage  SEM  work.  No  striking  difference  is  seen  between  cell  (a)  and  cells  (c)  and  (d).  In  the 
current  study  the  cell  performance  of  all  cells  is  better  than  the  performance  observed  previously 
(Technical  report  submitted  to  Optomec  Inc.  /AFRL,  October  2010).  However  the  impact  of  in¬ 
line  mixer  remains  unclear  and  needs  further  clarification/investigation. 
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Figure  36  Cell  (a)  -no  anode  interlayer,  (c)  and  (d)-cells  with  anode  interlayer 

4.2  Graded  Braze  Power  Electronics 

It  was  determined  after  showing  some  progress  with  silver  nano  paste  but  not  having  the  micro 
pen  operational  that  we  should  try  different  routes  to  achieve  two  tone  metallic  brazing.  We 
decided  to  pursue  a  tape  casting  route. 

Tape  caste  of  silver  nano  particles  was  performed  with  success  allowing  for  a  solid  tape  that 
could  be  handled,  hole  punched  and  used  as  the  center  of  the  die  attach  layer.  The  first 
experiments  used  4  single  layers  of  Ag  tape  as  the  die  attach  layer.  The  results  were  documented 
into  a  power  point  while  one  of  the  images  is  shown  in  Figure  39.  The  next  series  of  experiment 
used  the  hole-punched  Ag  tape  as  the  center  of  the  two  tone  layer.  Unfortunately  the  escaping 
gasses  skewed  the  layer  resulting  in  a  failed  die  attach.  It  has  been  determined  that  a  small  jig 
should  be  fabricated  that  will  constrain  the  sample  and  allow  for  an  even,  consistent,  and 
adequate  pressure  to  be  applied  during  the  sintering  process.  Along  with  developing  a  jig  it  was 
also  determined  that  a  material  set  be  purchased  that  would  mimic  the  actual  die  attach  scheme. 
Purchasing  of  DBC,  Sic  and  effective  back  side  metallization  of  the  Sic  wafers  was  investigated. 
Sources  have  been  found  and  a  proper  back  side  metallization  has  been  determined  (Ti,  Ni,  Au). 
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Figure  37  Cross  Section  of  a  Typical  Cell  with  Anode  Interlayer 


(a)  (c)  (d) 


Figure  38  Scanning  Electron  Image  of  Cells  (a),  (c),  and  (d) 
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Figure  39  As  Sintered  @320C  Ag  Tape  Die  Attach  Layer 


While  waiting  for  the  finalization  of  the  material  set,  side  experiments  have  been  investigated.  It 
was  mentioned  in  a  paper  that  finding  a  way  to  increase  the  oxygen  content  in  the  die  attach  layer 
may  help  in  the  decomposition  of  the  plastic  deep  in  the  center  of  the  layer.  It  was  thought  that 
adding  CuO  to  the  Ag  tape  may  supply  additional  oxygen  while  leaving  behind  copper  metal 
during  decomposition  at  elevated  temperatures.  A  tape  cast  of  CuO  was  made  and  fired  at  300  °C 
to  determine  if  the  plastic  decomposition  would  reduce  the  CuO  into  Cu  during  sintering.  This 
did  occur  but  the  material  quickly  returned  to  CuO  before  the  Cu  could  stabilize  due  to  being  in 
an  02  environment  at  elevated  temperatures.  From  this  result  it  was  detennined  that  even  if  the 
CuO  expedited  the  decomposition  of  the  plastic  it  would  re-oxidize  and  be  ineffective  as  a 
thermal  conductor  and  bonding  material. 

The  idea  of  the  addition  of  a  flux  has  been  tossed  around  so  incorporating  a  flux  into  the  tape 
was  investigated.  The  most  probable  chemical  reaction  that  could  occur  between  the  Ag  nano 
particles  with  the  plastic  tape  and  the  environment  were  determined  to  be  oxidation  of  the  Ag 
particles.  The  Gibbs  free  energy  of  this  reaction  was  looked  at  within  the  processing  temperature 
range.  It  was  detennined  that  after  approximately  200C  Ag  oxidation  ceases.  Knowing  this, 
future  investigation  into  the  synthesis  of  the  tape  caste  has  been  performed  looking  for  ways  to 
protect  the  Ag  nano  particle  from  oxidation  until  it  is  above  200  °C.  It  has  been  determined  that 
the  dispersant  will  play  the  biggest  role  in  preventing  oxidation.  This  has  led  to  looking  at  a  few 
different  tape  caste  recipes  that  use  a  more  aggressive  dispersant  with  strong  surface  absorption 
characteristics. 
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5  Measurement  of  Static  and  Dynamic  Performance 
Characteristics  of  Small  Electric  Propulsion  Systems 


Unmanned  aerial  vehicles  are  being  utilized  by  numerous  groups  around  the  world  for  various 
missions.  Most  of  the  smaller  vehicles  that  have  been  developed  use  commercially-off-the-shelf 
parts,  and  little  information  about  the  perfonnance  characteristics  of  the  propulsion  systems  is 
available  in  the  archival  literature.  In  light  of  this,  the  aim  of  the  present  research  was  to 
determine  the  performance  of  various  small-scale  propellers  in  the  4.0  to  6.0  inch  diameter  range 
driven  by  an  electric  motor.  An  experimental  test  stand  was  designed  and  constructed  in  which 
the  propeller/electric  motor  was  mounted  in  a  wind  tunnel  for  both  static  and  dynamic  testing, 
and  the  results  were  compared  to  those  from  previous  studies.  For  static  testing,  the  coefficient  of 
thrust,  the  coefficient  of  propeller  power,  and  the  total  propulsive  efficiency,  defined  as  the  ratio 
of  the  propeller  output  power  to  the  electrical  input  power,  were  plotted  versus  the  propeller 
rotational  speed.  For  dynamic  testing,  the  rotational  speed  of  the  propeller  was  held  constant  at 
regular  intervals  while  the  airspeed  was  increased  from  zero  to  the  windmill  state.  The 
coefficient  of  thrust,  the  coefficient  of  propeller  power  and  the  propeller  efficiency  were  plotted 
versus  the  advance  ratio  for  various  rotational  speeds.  The  thrust  and  torque  were  found  to 
increase  with  rotational  speed,  propeller  pitch  and  diameter,  and  decrease  with  airspeed.  Using 
the  present  results  and  data  from  archival  and  non-archival  sources,  it  was  found  that  the 
coefficient  of  thrust  could  not  be  correlated  with  propeller  diameter  for  square  propellers  where 
D  =  P.  For  a  family  of  propellers  (same  manufacturer  and  application),  correlations  for  the 
coefficient  of  thrust,  the  coefficient  of  propeller  power  and  the  propeller  efficiency  could  be 
improved  by  modifying  either  the  coefficients  or  the  advance  ratio  with  DIP.  This  dimensionless 
ratio  allows  for  the  propeller  pitch  to  be  accounted  for  in  the  perfonnance  coefficients. 

5.1  Introduction 

Interest  in  the  performance  of  small  propellers  operating  at  low  Reynolds  numbers  has  grown 
recently.  The  aerospace  industry  has  developed  numerous  unmanned  aerial  vehicles  (UAVs)  and 
has  kept  most  of  the  data  about  the  propulsion  systems  proprietary.  Very  little  information  is 
available  in  the  archival  literature  about  the  performance  characteristics  of  these  motor  and 
propeller  combinations.  The  present  research  and  others  like  it  have  aimed  to  gather  and  compare 
information  about  these  small  propulsion  systems  so  that  proper  motor  and  propeller 
combinations  can  be  selected  for  a  given  mission  profile.  Several  papers  were  reviewed  that 
relate  directly  to  the  present  work  and  provide  direction  for  the  research. 

Brandt  and  Selig  experimentally  determined  efficiency  as  well  as  coefficients  of  thrust  and 
power  for  low  Reynolds  number  propellers  [85].  The  parametric  ranges  were  as  follows: 

Propeller  diameter  9  <  D  <  1 1  inches,  propeller  rotational  speed  1500  <  n  <  7500  RPM,  and  the 
incoming  air  velocity  V<^  ranged  from  zero  (static)  to  the  windmill  state  of  each  propeller,  i.e., 
that  point  at  which  the  propeller  generates  zero  thrust.  A  test  stand  was  built  inside  the  UIUC 
wind  tunnel  to  measure  thrust,  torque,  and  propeller  rotational  speed.  Freestream  air  velocity  was 
measured  using  a  Pitot  tube  and  one  of  two  differential  pressure  transducers  depending  on  the 
airspeed  range.  Velocity  corrections  were  applied  to  account  for  the  change  in  upstream  airspeed 
at  the  Pitot  tube  created  by  the  propeller  as  well  as  the  pressure  change  created  by  the  fairing  and 
the  constriction  of  the  propeller  slipstream  caused  by  the  walls.  In  total,  79  propellers  from  four 
different  manufacturers  were  tested  to  find  the  coefficient  of  thrust,  the  coefficient  of  power  and 
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the  propeller  efficiency,  all  of  which  were  plotted  against  advance  ratio.  The  designs  of  the 
propellers  ranged  from  those  for  electric  motors  to  those  used  for  fuel-powered  engines.  For  each 
test,  the  rotational  speed  of  each  propeller  was  fixed  while  the  freestream  airspeed  was  varied. 
Four  different  values  of  propeller  rotational  speed  (n  =  3000,  4000,  5000,  and  6000  RPM)  were 
tested  for  each  of  the  propellers.  The  results  show  that  the  propeller  efficiency  increases  with  the 
propeller  speed.  This  is  primarily  due  to  the  increase  in  Reynolds  number  as  the  propeller  spins 
faster.  Overall,  the  propeller  efficiency  ranged  from  28  <  //p  <  65%.  The  propellers  were  also 
tested  statically,  but  the  data  is  only  available  in  the  UIUC  propeller  database  [96]. 

Gamble  designed  an  intricate  Lab  VIEW  program  to  automatically  collect  data  and  generate 
propeller  perfonnance  plots  [89].  A  dynamometer  was  constructed  using  beam-type  load  cells  to 
measure  thrust  and  torque.  The  development  of  the  Lab  VIEW  program  was  detailed  as  well  as  a 
procedure  for  carrying  out  the  experiment.  Propellers  were  tested  for  repeatability  by  perfonning 
identical  experiments  over  several  days  with  two  identical  propellers.  The  results  primarily  focus 
on  the  effect  of  the  Reynolds  number  on  thrust  and  power  coefficients  and  efficiency  versus 
advance  ratio.  Thrust  versus  velocity  was  compared  for  propellers  with  constant  diameter  and 
varying  pitch.  Lastly,  advance  ratio  was  modified  by  replacing  diameter  with  pitch  in  the 
equation  for  advance  ratio.  The  optimal  advance  ratio  is  shown  using  this  technique.  This  allows 
for  the  optimal  pitch  of  a  model  propeller  to  be  selected  to  achieve  maximum  efficiency.  The 
diameter  can  then  be  chosen  from  plots  of  thrust  versus  velocity  to  produce  the  required  thrust 
for  the  airframe. 

Deters  and  Selig  perfonned  static  tests  on  smaller  propellers  ranging  from  2.5  <  D  <  5  inches  in 
diameter  [88].  Static  coefficients  of  thrust  and  power  as  well  as  the  figure  of  merit  (FOM  = 
C^2/V2CP,  typically  used  to  measure  the  efficiency  of  helicopters)  using  modified  coefficients 
of  thrust  and  power  that  use  disk  area  and  tip  speed  were  detennined  experimentally.  The  test 
stand  utilized  a  0.3  kg  load  cell  and  a  25  oz-in  torque  transducer  to  measure  thrust  and  torque, 
respectively.  Propeller  rotational  speeds  ranging  from  2500  <  n  <  27,000  RPM  were  measured 
using  an  infrared  detector.  A  schematic  of  the  test  stand  indicated  the  locations  of  the 
components  and  a  fairing  surrounding  the  load  cell  and  torque  transducer.  Calibrations  of  the 
components  were  performed  and  data  was  collected  using  a  data  acquisition  board.  The  geometry 
of  each  propeller  was  found  using  PropellerScanner  software  to  find  the  chord  and  twist 
distribution  [92],  This  was  used  to  calculate  the  Reynolds  number  at  the  75%  chord  location. 
Results  show  that  over  the  rotational  speed  range  tested,  the  figure  of  merit  remained  fairly 
constant  throughout  the  test.  The  results  also  show  that  a  larger  diameter  propeller  is  more 
efficient  than  a  smaller  one,  and  a  propeller  with  a  lower  pitch  is  more  efficient  than  one  with  a 
higher  pitch. 

01  et  al.  took  a  more  analytical  approach  to  studying  small  propellers  operating  at  low  Reynolds 
numbers  [95].  Iterative  methods  were  used  to  calculate  the  coefficient  of  thrust,  the  coefficient  of 
torque,  and  the  propeller  efficiency  using  propeller  momentum  theory  and  blade-element 
methods.  Propellers  were  discretized  by  cutting  and  tracing  sections  as  well  as  digital  scans. 
Leading  and  trailing  edges  were  fitted  to  the  UIUC  propeller  library  so  that  the  resulting  analysis 
in  XFOIL  would  successfully  converge.  The  iterative  process  for  thrust  was  dependent  on  the 
various  Reynolds  numbers  across  the  propeller  blade  at  a  given  rotational  speed.  Two  separate 
experimental  setups  were  constructed  to  compare  the  numerical  results.  Propellers  in  the  6  <D< 
12  inch  range  were  tested  in  the  Langley  Research  Center  Basic  Aerodynamics  Research  Tunnel 
(BART)  and  larger  propellers  in  the  14  <  D  <  20  inch  range  were  tested  in  the  AFRL  Vertical 
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Wind  Tunnel  (VWT).  Static  tests  were  performed  with  the  wind  tunnel  sides  open  to  alleviate  the 
induced  airflow  velocity  inside  the  wind  tunnel.  Blockage  corrections  were  applied  to  BART 
tests  but  not  to  VWT  tests,  since  the  tunnel  diameter  of  the  VWT  was  greater  than  five  times  the 
diameter  of  the  propellers  tested.  Drag  on  the  test  stand  was  corrected  by  sweeping  tunnel 
velocity  and  generating  curve  fits  that  were  used  to  adjust  the  actual  data.  A  large  sensitivity  to 
twist  distribution  was  observed  in  the  tests  and  the  analysis.  01  et  al.  postulated  that  plots  of 
torque  coefficient  versus  advance  ratio  are  sometimes  misleading  because  they  do  not  account 
for  Reynolds  number  effects.  It  was  also  shown  that  when  the  ratio  of  diameter  to  pitch  is  scaled 
(10  x  10  to  12  x  12,  for  example)  the  experimental  data  fits  together  well  within  the  bounds  of 
error.  Modifications  to  the  dimensionless  tenns  to  factor  in  propeller  pitch  were  presented, 
however  more  research  was  deemed  necessary  to  apply  this  theory. 

Corrigan  and  Altman  examined  different  methods  for  wind  tunnel  blockage  corrections  [87]. 
These  methods  included  the  Glauert  correction  as  well  as  a  correction  by  Hackett  et  al.  [90,  91]. 
These  methods  were  described  in  detail  and  their  applications  were  shown.  A  wind  tunnel 
experiment  was  designed  and  constructed  to  record  the  necessary  variables  to  calculate  total 
propulsive  efficiency.  This  is  in  contrast  to  other  works  that  primarily  explored  propeller 
efficiency.  The  stand  was  constructed  using  a  beam-type  load  cell  and  a  reaction  torque  sensor. 
Three  propellers  ( D  =  10,  12,  and  14  inches)  were  tested  using  different  motors  for  each 
propeller.  Static  pressure  taps  were  used  on  the  wall  of  the  wind  tunnel  test  section  to  record  the 
changes  in  pressure  forward  and  aft  of  the  propeller  disk  plane  for  the  velocity  corrections.  The 
Glauert  method  did  not  provide  sufficient  correction  for  large  blockage  conditions.  The  Hackett 
method  yielded  more  correction  at  higher  airspeeds  and  larger  propeller  diameters,  but  the 
method  could  not  be  validated  and  therefore  further  work  was  found  to  be  necessary. 

Merchant  and  Miller  perfonned  dynamic  tests  on  propellers  in  the  6  <  D  <  22  inch  range  [94],  A 
test  stand  was  constructed  to  record  propeller  performance  parameters,  where  the  thrust  and 
torque  were  collected  by  a  combined  thrust/torque  cell.  The  load/torque  cell  was  calibrated  using 
dead  weights  in  the  axial  (thrust)  and  transverse  (torque)  directions.  Wind  tunnel  velocity  was 
measured  directly  using  a  Pitot  probe  and  a  differential  pressure  transducer.  Since  the  propellers 
were  large  compared  to  the  test  section,  blockage  corrections  developed  by  Glauert  were  applied 
to  the  results.  Readings  were  taken  at  wind-off-zero  conditions  before  and  after  each  test  [90]. 
These  values  were  then  averaged  and  subtracted  from  the  test  data  to  account  for  zero  drift  and 
temperature  effects.  Data  was  collected  at  constant  propeller  rotational  speeds  and  the  wind 
tunnel  velocity  was  varied  to  sweep  through  values  of  advance  ratio.  The  results  were  compared 
to  other  works  and  were  shown  to  be  acceptable.  The  setup  was  also  tested  for  variations  in  flow 
angularity.  Pitch  and  yaw  variations  between  -3  and  +3  arc  degrees  were  examined  and  it  was 
shown  that  only  the  coefficient  of  thrust  was  affected  by  a  change  in  pitch.  However,  it  was 
shown  that  pitch  variations  of -3  and  +3  degrees  yielded  the  same  results,  which  indicated  that 
the  system  was  symmetric  in  the  pitch  direction.  Lastly,  two  identical  propellers  made  by  the 
same  manufacturer  were  tested  and  compared,  which  showed  that  for  some  propellers  there  may 
be  significant  differences  in  performance  due  to  manufacturing.  Very  limited  results  were 
presented,  however,  and  the  results  shown  only  give  a  small  sample  of  the  entire  test  range. 

The  objective  of  the  present  research  was  to  determine  the  performance  of  various  commercially- 
available  small-scale  propellers  driven  by  an  electric  motor.  An  experimental  test  stand  was 
designed  and  constructed  in  which  the  electric  motor  was  mounted  in  a  wind  tunnel  at  Wright 
State  University  for  both  static  and  dynamic  testing.  The  freestream  airspeed  was  varied  from 
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zero  to  the  windmill  state  for  each  propeller.  The  rotational  speed  was  varied  over  the  operational 
range  recommended  by  the  propeller  manufacturers,  while  ensuring  that  the  electric  motor  did 
not  overheat.  The  primary  measurement  devices  were  calibrated,  and  an  extensive  uncertainty 
analysis  was  perfonned.  The  results  from  the  present  experiment  were  compared  to  those  from 
previous  studies  for  both  static  and  dynamic  data.  For  static  testing,  the  coefficient  of  thrust,  the 
coefficient  of  propeller  power  and  the  total  propulsive  efficiency  were  plotted  versus  the 
propeller  rotational  speed.  For  dynamic  testing,  the  rotational  speed  of  the  propeller  was  held 
constant  at  regular  intervals  while  the  freestream  airspeed  was  increased  from  zero  to  the 
maximum.  The  coefficient  of  thrust,  the  coefficient  of  propeller  power  and  the  propeller 
efficiency  were  plotted  versus  the  advance  ratio  for  various  rotational  speeds. 

5.2  Background 

The  performance  characteristics  to  be  detennined  by  the  experimental  setup  are  as  follows.  The 
coefficients  of  thrust,  torque  and  propeller  power,  and  the  propeller  efficiency  are  Merchant  and 
Miller  [94]: 
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The  three  perfonnance  coefficients  and  the  propeller  efficiency  defined  above  are  typically 
plotted  against  the  advance  ratio  for  dynamic  testing: 


(5-1) 


The  Glauert  correction  variable  is: 
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The  propeller  disk  area  and  wind  tunnel  area  are,  respectively: 

Tr  D  ^ 

AP  =  ^~,  Awt  =  W2  (5-5) 

The  corrected  thrust  is  defined  as  the  measured  thrust  minus  the  drag  force  due  to  the  flow  of  air 
over  the  motor,  torque  cell  and  load  cell  [97]: 

T'  =  T-Fd  (5-6) 
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The  total  propulsive  efficiency  is  the  ratio  of  the  propeller  output  power  to  the  electrical  input 
power: 
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(5-7) 


The  density  of  air  is  given  by  the  perfect  gas  law: 
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(5-8) 


5.3  Experimental  setup 

The  objective  of  the  present  experiment  was  to  determine  the  performance  characteristics  of 
small  electric  motor/propeller  combinations  from  static  conditions  to  the  windmill  state.  The 
overall  design  of  the  dynamic  test  rig  is  shown  in  Figure  40.  The  electric  motor  was  directly 
attached  to  a  25  oz-in  torque  cell  (Transducer  Techniques,  Model  RTS-25),  which  was  able  to 
withstand  10  kg  in  thrust  and  1.7  kg  in  shear.  The  torque  cell  was  in  turn  mounted  onto  a  1-kg 
single  point  beam-type  load  cell  (Transducer  Techniques,  Model  LSP-1).  Each  cell  was  driven 
by  a  signal  conditioner  (Transducer  Techniques,  Model  TMO-1)  that  produced  a  0  to  5  Volt 
linear  output.  The  assembly  of  the  motor,  torque  cell  and  load  cell  is  shown  in  Figure  41.  The 
motor  was  held  in  place  with  a  custom-designed  clam-shell  clamp,  in  which  fins  were 
incorporated  to  increase  the  convective  heat  transfer  from  the  electric  motor. 


Figure  40  Schematic  Diagram  of  the  Experimental  Setup 
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(a) 


(b) 


Figure  41  Assembly  of  Motor,  Torque  Cell  and  Load  Cell:  (a)  Solid  Model  Representation;  (b) 

Photograph 

The  load  cell  was  attached  to  a  section  of  1.25-inch  square  aluminum  tubing,  which  acted  as  a 
riser  to  place  the  propeller  in  the  middle  of  the  test  section.  The  bottom  of  the  riser  was 
connected  to  an  optical  breadboard  table  (Melles-Griot,  Model  BBSS-25-610-1219)  using 
flanges  of  angle  aluminum.  A  hole  was  milled  in  the  acrylic  floor  of  the  wind  tunnel  for  the 
aluminum  riser  to  pass  through.  The  low-speed  wind  tunnel  at  Wright  State  University  is  an  open 
circuit  design  capable  of  producing  speeds  from  0.6  to  36  m/s  with  a  contraction  ratio  of  6.25:1. 
The  square  entrance  of  the  wind  tunnel  has  a  3.8  m2  opening  with  aluminum  hexagonal 
honeycomb  sections  that  serve  as  a  flow  straightener.  The  height  and  width  of  the  square  test 
section  is  W=  0.6096  m,  and  its  length  is  2.438  m.  Doors  on  one  side  of  the  test  section  allow  for 
an  entire  wall  to  be  opened  for  easy  access. 

The  data  acquisition  system  used  to  collect  data  from  the  instrumentation  consisted  of  a  DAQ 
board  (National  Instruments,  Model  SCC-68)  and  a  DAQ  card  (National  Instruments,  Model 
PCI-6221)  installed  in  a  PC.  Shielded  wires  were  used  to  connect  the  outputs  of  the  transducers 
to  the  DAQ  board.  The  electric  motor  driving  the  propeller  was  energized  using  a  precision  DC 
power  supply  (Hewlett-Packard,  Model  6012B).  A  servo  tester  (GWS,  Model  MT-1)  was  used  to 
control  the  rotational  speed  of  the  propeller.  The  voltage  supplied  to  the  electric  motor  was 
measured  using  a  digital  multi-meter  (National  Instruments,  Model  USB-4065).  To  measure  the 
current,  a  DC  Hall  effect  current  transducer  (CR  Magnetics,  Model  CR52 10-30)  with  a  range  of 
0  to  30  A  was  placed  in-line  between  the  power  supply  and  the  motor  speed  controller. 

A  remote  optical  sensor  (Monarch  Instrument,  Model  ROS-W)  connected  to  a  panel  meter 
(Monarch  Instrument,  Model  ACT-3X)  was  used  to  measure  propeller  rotational  speed. 
Reflective  tape  supplied  with  the  sensor  was  placed  near  the  hub  on  the  leeward  side  of  the 
propeller  so  that  the  optical  sensor  did  not  have  to  be  adjusted  between  runs. 

Atmospheric  pressure  was  measured  to  determine  the  density  of  the  air.  To  record  atmospheric 
pressure,  a  barometer  (Vaisala,  Model  PTB1 10)  capable  of  measuring  500  to  1100  mbar  with 
accuracy  of  ±0.3  mbar  was  used.  The  differential  pressure  produced  by  the  Pitot  tube  was 
measured  using  a  differential  pressure  manometer  (MKS,  Model  226A).  The  height  of  the  Pitot 
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tube  from  the  floor  of  the  wind  tunnel  was  selected  by  traversing  the  boundary  layer  thickness 
using  the  Pitot  tube  as  outlined.  The  height  was  set  to  H  =  2.5  inches,  and  the  Pitot  tube  was 
made  parallel  to  the  wind  tunnel  walls  by  using  a  bubble  level  and  a  custom-made  jig. 

The  temperature  of  the  motor  was  measured  using  a  Type  T  thermocouple  while  the  temperature 
of  the  air  inside  the  wind  tunnel  was  measured  using  a  Type  E  thermocouple  probe  (Omega, 
Model  EMQSS-125G-12).  The  Type  T  thermocouple  junction  was  placed  on  the  center  of  the 
motor  and  was  held  in  place  by  the  aluminum  clam-shell  clamp.  The  Type  E  probe  was  mounted 
in  the  floor  of  the  wind  tunnel  ahead  of  the  motor/propeller  so  that  the  sensing  junction  extended 
into  the  airflow.  The  thermocouples  were  connected  to  thermocouple  modules  (National 
Instruments,  Model  SCC-TC01)  on  the  data  acquisition  board.  The  signals  from  the  eight  sensors 
were  read  using  custom-designed  Lab  VIEW  virtual  instruments. 

The  twenty-four  propellers  selected  for  analysis  ranged  from  4.0  <  D  <  6.0  inches  in  diameter 
and  2.0  <  P  <  5.5  inches  in  pitch,  as  shown  in  Table  7.  Some  of  the  propellers  were  selected  to 
overlap  with  previous  research  so  that  the  procedures  and  test  setup  used  for  the  measurements 
could  be  compared  and  validated.  The  GWS  4.5  x  3.0  and  GWS  5.0  x  4.3  inch  propellers  were 
tested  statically  and  compared  to  Deters  and  Selig  [88],  An  APC  8.0  x  3.8  inch  Slow  Flyer  was 
tested  dynamically  and  compared  to  the  results  posted  on  the  UIUC  Propeller  Database  while  an 
APC  6.0  x  4.0  inch  propeller  was  tested  dynamically  and  compared  to  the  results  presented  by  01 
et  al.  [95,  96]. 
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Table  7  Summary  of  Propellers  Studied. 


Manufacturer 

Nominal  D/P  (in  x  int 

Desipnation 

APC 

4.10  x  4.10 

Sneed  400  Electric 

APC 

4.20  x  2.00 

Snort 

APC 

4.20  x  4.00 

Free  Flight 

APC 

4.50  x  4.10 

Sneed  400  Electric 

APC 

4.70  x  4.25 

Sneed  400  Electric 

APC 

4.75  x  4.75 

Sneed  400  Electric 

APC 

4.75  x  5.50 

Sneed  400  Electric 

APC 

5.10  x  4.50E 

Thin  Electric 

APC 

5.25  x  4.75 

Sneed  400  Electric 

APC 

5.50  x  2.00 

Free  Flight 

APC 

5.50  x  4.50 

Sneed  400  Electric 

APC 

6.00  x  2.00 

Snort 

APC 

6.00  x  4.00  E 

Sneed  400  Electric 

APC 

8.00  x  3.8 

Slow  Fiver 

Graunner 

4.00  x  3.00 

Cam  Sneed 

Graunner 

4.70  x  4.00 

Cam  Sneed 

Graunner 

4.70  x  4.70 

Cam  Sneed 

Graunner 

5.50  x  4.30 

Cam  Sneed 

Graunner 

5.50  x  5.50 

Cam  Sneed 

GWS 

4.00  x  2.50 

GWS 

4.00  x  4.00 

GWS 

4.50  x  3.00 

GWS 

5.00  x  3.00 

GWS 

5.00  x  4.30 

5.4  Uncertainty  analysis 

The  uncertainties  of  all  of  the  calculated  results  described  in  the  above  equations  were 
detennined  using  the  root-sum-square  uncertainty  method  [93].  During  experimentation,  eight 
primary  measurements  were  made  using  the  data  acquisition  system:  Uncorrected  thrust  (T); 
torque  (Q):  propeller  rotational  speed  (n);  atmospheric  pressure  (Patm);  atmospheric  temperature 
(f at m) ;  Pitot  tube  pressure  difference  (PdlfT);  motor  voltage  (V);  and  motor  amperage  (/).  The 
uncertainty  of  a  given  measurement  was  estimated  to  be  the  sum  of  the  calibration  uncertainty 
and  the  confidence  interval  of  the  collected  data  set  at  a  confidence  level  of  99%: 

AU  =  AUcal  +  A  U99  (5-9) 
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In  order  to  calibrate  the  torque  cell,  two  identical  arms  were  attached  to  the  sides  of  the  motor 
clamp  so  that  the  torque  cell  could  be  calibrated  in  both  directions  of  rotation  simultaneously. 
Varying  weights  were  hung  from  one  of  the  arms  to  calibrate  in  the  clockwise  direction,  and  then 
the  process  was  repeated  for  the  counterclockwise  direction.  The  load  cell  used  to  measure  thrust 
was  calibrated  in  situ  as  follows:  A  strand  of  fishing  line  was  attached  to  the  front  of  the 
propeller  using  aircraft  wire.  This  strand  was  then  passed  over  a  smooth  cylinder  with  bearings 
mounted  in  the  wind  tunnel.  Varying  weights  were  suspended  from  the  fishing  line  over  the 
expected  range  of  thrust,  and  voltage  readings  were  recorded. 

The  drag  of  the  fixture  was  measured  versus  airspeed  by  removing  the  propeller  and  replacing  it 
with  a  propeller  hub  with  the  blades  removed.  The  airspeed  was  increased  systematically  while 
data  was  collected  from  the  load  cell  and  the  Pitot  tube.  The  free-stream  velocity  was  then 
calculated  and  the  measured  drag  was  plotted  against  the  velocity.  A  second-order  regression 
was  applied  to  the  points  and  this  equation  was  used  in  the  calculation  of  the  corrected  thrust. 

Table  8  gives  the  uncertainties  for  each  device  or  transducer  used  to  collect  the  data.  The 
principal  equations  used  for  detennining  the  uncertainties  of  the  computed  quantities  shown  in 
the  graphs  in  the  Results  and  Discussion  section  are  shown  below. 

Coefficient  of  Thrust: 


(5-10) 


Coefficient  of  Torque: 


A  CQ  = 


+ 


1 


2-12 


(5-11) 


Coefficient  of  Power: 


A  CP  = 


+ 


1 

2 


(5-12) 
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Propeller  Efficiency. 


Ary?  = 


1 

l~JCjACPJ* 


(5-13) 


Advance  Ratio : 


A /  = 


MC\2  /-CAn\2  /-CAD\2!^ 

\nD )  n2D  )  nD2  ) 


(5-14) 


Total  Propulsive  Efficiency. 
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Table  8  Uncertainties  of  Primary  Measurement  Sensors  and  Calibration  Sources. 


Measurement 

Sensor 

Uncertainty 

Thrust,  T 

Transducer  Techniques  LSP  1kg  Load 
Cell 

AFcai  =  ±7.70  g 

Torque,  Q 

Transducer  Techniques  RTS  25  oz-in 
Reaction  Torque  Sensor 

A<9cai  =  ±0.0498  g- 
m 

Atmospheric  Temperature, 

Fatm 

Omega  Type  E  Thermocouple 

AFatai.cal  — 

±0.0334  °C 

Calibration  Mass 

Ohaus  Digital  Scale 

A m  =  ±1.00  x  10"3 
g 

Propeller  Diameter,  D 

Digital  Vernier  Calipers 

AD  =  ±1.00  x  10'5 
m 

Propeller  Rotational  Speed,  n 

Monarch  Instruments  Remote  Optical 
Sensor  (ROS)  and  ACT  3x  Panel 
Tachometer 

An  =  ±  1  RPM 

Motor  Voltage,  V 

National  Instruments  USB-4065  Digital 
Multi-Meter 

AV=±  1.00  x  10'3 

V 

Motor  Current,  / 

CR  Magnetics  CR52 10-30  Current 
Transducer 

A/=  ±  (1%  x 
Reading) 

Atmospheric  Pressure,  Pabs 

Vaisala  PTB1 10  Barometer 

APatm  =  ±  30.0  Pa 

Pitot  Tube  Differential 
Pressure,  /fim- 

MKS  226A  Differential  Pressure 
Manometer 

A  Pd  iff  =  ±  (0.3%  x 
Reading) 

5.5  Experimental  procedures 

Two  separate  procedures  were  developed  for  the  static  and  dynamic  tests.  For  all  of  the  tests,  the 
power  supply  driving  the  motor  controller  for  the  propeller  motor  was  turned  on  and  set  to  a 
nominal  output  of  1 1.1  Volts,  which  matches  the  voltage  output  of  a  standard  3-cell  lithium- 
polymer  battery.  Then,  the  data  acquisition  system  and  the  signal  conditioners  driving  the 
sensors  were  powered  up  for  the  warm-up  periods  recommended  by  the  manufacturers. 

5.5.1  Static  test  procedure 

After  the  warm-up  period,  the  load  cell  and  torque  cell  were  zeroed  by  adjusting  the  balance 
potentiometers  on  the  signal  conditioners  so  that  the  voltage  outputs  were  as  close  as  possible  to 
zero.  At  this  point,  five  hundred  data  points  were  collected  with  the  propeller  motor  off  to  obtain 
baseline  values  for  the  load  cell  and  torque  cell.  The  propeller  motor  was  then  set  to  the  first 
desired  speed  setting  and  one  thousand  data  points  were  collected.  The  propeller  motor  was  then 
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turned  off  and  another  set  of  500  data  points  was  acquired.  The  average  values  for  thrust  and 
torque  from  the  two  propeller-off  states  were  averaged  and  this  value  was  used  to  correct  the 
thrust  and  torque  measurements  to  account  for  zero  drift  and  temperature  effects  (Merchant  and 
Miller,  2008).  The  process  was  then  repeated  for  increased  values  of  rotational  speed  until  the 
maximum  speed  was  achieved. 

5.5.2  Dynamic  test  procedure 

After  the  wann-up  period,  the  differential  pressure  transducer  reading  the  Pitot  tube  and  the 
signal  conditioners  reading  the  load  cell  and  the  torque  cell  were  zeroed.  Five  hundred  data 
points  were  taken  with  the  propeller  motor  off  and  the  wind  tunnel  motor  off.  At  the  end  of  the 
first  five  hundred  points,  the  propeller  motor  was  set  to  the  desired  rotational  speed  setting  and 
the  wind  tunnel  airspeed  was  set  to  the  first  desired  setting.  After  the  system  reached  steady  state, 
one  thousand  data  points  were  acquired.  Next,  the  wind  tunnel  airspeed  setting  was  changed  and 
the  propeller  rotational  speed  was  adjusted  to  match  the  original  setting.  This  process  was 
repeated  until  the  windmill  state  of  the  propeller  was  reached.  The  propeller  motor  and  the  wind 
tunnel  motor  were  both  stopped  at  this  point,  and  then  five  hundred  data  points  were  collected  in 
order  to  again  account  for  drift  in  the  sensors.  Data  sets  were  collected  for  approximately  ten 
wind  tunnel  airspeed  settings  for  each  of  the  four  rotational  speed  settings  for  each  propeller 
tested. 

5.6  Results  and  discussion 

To  ensure  that  the  collected  data  was  repeatable  and  correct,  tests  were  necessary  to  validate  the 
static  and  dynamic  results.  The  first  type  of  test  checked  for  repeatability  of  the  same  propeller  as 
well  as  the  repeatability  across  three  identical  propellers.  The  second  type  of  test  was  to  compare 
the  results  of  the  present  experiment  to  published  results  from  researchers  using  the  same 
propeller.  A  complete  summary  of  the  data  collected  from  the  static  and  dynamic  tests  is 
provided  by  Brezina  [86]. 

5.6. 1  Validation  of  the  static  test 

To  check  the  repeatability  of  the  experiment,  three  identical  Graupner  4.7  x  4.7  inch  propellers 
were  tested  under  static  conditions  three  times  each,  thus  creating  a  total  of  nine  sets  of  data. 

This  was  done  to  determine  the  repeatability  of  the  experiment  for  multiple  runs  of  the  same 
propeller  as  well  as  establishing  whether  manufacturing  variability  affected  the  performance  of 
identical  propellers.  Figure  42  shows  typical  results  for  a  static  propeller,  where  both  the  thrust 
and  torque  increase  monotonically  with  rotational  speed.  The  coefficients  of  thrust  and  propeller 
power  were  relatively  constant,  whereas  the  total  efficiency  reached  a  peak  value  at 
approximately  n  =  20,000  RPM.  The  uncertainties  of  the  coefficients  of  thrust  and  propeller 
power  increased  significantly  at  the  lowest  propeller  rotational  speed.  This  was  driven  by  the 
uncertainty  of  the  load  cell  and  the  torque  cell  at  relatively  small  values  of  thrust  and  torque. 
Figure  43  shows  that  the  repeatability  of  the  reduced  data  (coefficient  of  thrust,  coefficient  of 
propeller  power  and  total  propulsive  efficiency)  was  excellent.  The  data  from  all  nine  tests  fall 
within  the  uncertainty  bounds  for  the  first  run.  The  duplicate  propellers  also  fall  directly  in  line, 
meaning  that,  at  least  for  this  type  of  propeller,  manufacturing  differences  can  be  neglected. 
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Figure  42  Typical  Static  Test  Results 

(Graupner  4.7  x 4.7  inch  Propeller):  (a)  Thrust  and  Torque  versus  Rotational  Speed,  (b)  Coefficient  of 
Thrust,  Coefficient  of  Propeller  Power  and  Total  Propulsive  Efficiency  versus  Rotational  Speed. 
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a) 


(b) 


(c) 


Figure  43  Comparison  of  Three  Identical  Propellers  under  Static  Testing 
(Graupner  4.7  x  4.7):  (a)  Coefficient  of  Thrust,  (b)  Coefficient  of  Propeller  Power,  (c)  Total  Propulsive 

Efficiency. 
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Static  tests  were  performed  on  two  propellers  (GWS  4.5  x  3.0  and  GWS  5.0  x  4.3)  which 
matched  tests  performed  by  Deters  and  Selig  [88].  The  coefficient  of  thrust  and  the  coefficient  of 
power  were  compared  to  data  provided  by  Deters  and  Selig  as  shown  in  Figure  44,  where  the 
results  for  both  propellers  show  good  agreement. 
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Figure  44  Comparison  of  the  Present  Static  Results  to  Deters  and  Selig  (2008) 

( GWS  4.5  x  3.0  and  GWS  5.0  x  4.3  Propellers):  (a)  Coefficient  of  Thrust,  (b)  Coefficient  of  Propeller 

Power. 

5.6.2  Static  test  results 

Having  established  the  validity  of  the  experimental  results,  static  test  data  was  collected  for  all  of 
the  propellers  shown  in  Table  7.  Figure  45  shows  a  comparison  between  propellers  with  constant 
diameter  and  varying  pitch,  while  Figure  46  gives  a  comparison  between  propellers  with  varying 
diameter  and  constant  pitch.  In  Figure  45,  the  coefficients  of  thrust  and  propeller  power  were 
relatively  constant  while  the  propeller  efficiency  increased  with  propeller  rotational  speed. 
Increasing  the  pitch  (while  holding  the  propeller  diameter  constant)  significantly  increased  all 
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three  measures  of  performance.  This  same  trend  can  be  found  in  the  data  provided  by  Deters  and 
Selig  [88]  for  the  coefficient  of  thrust  and  coefficient  of  propeller  power  for  the  GWS  4.0  x  4.0 
propeller  versus  that  for  the  GWS  4.0  x  2.5  propeller.  In  Figure  46,  the  variation  of  the  three 
perfonnance  parameters  with  propeller  diameter  is  also  shown  to  be  significant,  where  increasing 
the  diameter  decreased  the  thrust  coefficient  and  the  propeller  power  coefficient  but  increased 
the  total  propulsive  efficiency.  This  trend  is  also  apparent  in  the  data  reported  by  Deters  and 
Selig  for  the  following  propellers:  GWS  3.0  x  3.0,  GWS  4.5  x  3.0,  and  GWS  5.0  x  3.0. 
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Figure  45  The  Effect  of  Varying  Pitch  While  Holding  Diameter  Constant  (Static  Testing) 
(a)  Coefficient  of  Thrust,  (b)  Coefficient  of  Propeller  Power,  (c)  Total  Propulsive  Efficiency. 
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Figure  46  The  Effect  of  Varying  Diameter  While  Holding  Pitch  Constant  (Static  Testing) 
(a)  Coefficient  of  Thrust,  (b)  Coefficient  of  Propeller  Power,  (c)  Total  Propulsive  Efficiency. 
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5.6.3  Validation  of  the  dynamic  test 

The  dynamic  test  procedure  and  experiment  were  validated  by  comparing  the  present 
results  to  non-archival  and  archival  literature  sources.  Figure  47  shows  typical  dynamic  results 
for  the  thrust  and  torque  generated  by  one  propeller  over  the  full  range  of  airspeed  and  various 
levels  of  rotational  speed.  Both  the  thrust  and  torque  increased  with  rotational  speed  and 
decreased  with  airspeed,  as  expected.  For  a  given  rotational  speed,  the  coefficients  of  thrust  and 
propeller  power  decreased  with  advance  ratio,  whereas  the  propeller  efficiency  increased  with 
advance  ratio.  The  APC  8.0  x  3.8  inch  Slow  Flyer  propeller  was  tested  at  nominal  propeller 
rotational  speeds  of  n  =  4000  and  7000  rpm,  and  the  results  for  the  coefficient  of  thrust,  the 
coefficient  of  propeller  power,  and  the  propeller  efficiency  versus  advance  ratio  were  compared 
to  those  reported  on  the  UIUC  propeller  database  [96],  as  shown  in  Figure  48.  The  agreement 
with  the  data  from  the  UIUC  database  is  excellent  for  both  rotational  speeds,  even  where  the 
propeller  efficiency  drops  off  steeply  with  advance  ratio. 
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Figure  47  Typical  Dynamic  Test  Results  (Graupner  4.7  x  4.7  inch  Propeller) 

(a)  Thrust  and  Torque  versus  Airspeed  for  Various  Rotational  Speeds,  (b)  Coefficient  of  Thrust, 
Coefficient  of  Propeller  Power  and  Propeller  Efficiency  versus  Advance  Ratio. 
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Figure  48  Comparison  of  Present  Dynamic  Results  to  Selig  (2012)  (APC  8.0  x  3.8  SF) 

(a)  Coefficient  of  Thrust,  (b)  Coefficient  of  Propeller  Power,  (c)  Propeller  Efficiency. 

To  further  validate  the  dynamic  results,  an  APC  6.0  x  4.0  inch  propeller  was  tested  at  nominal 
propeller  rotational  speeds  of  n  =  8000  to  16000  rptn  by  intervals  of  2000  rptn  and  compared  to 
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the  results  reported  by  01  et  al.  [95],  as  shown  in  Figure  49.  The  coefficients  of  thrust  and  torque 
decrease  with  advance  ratio  and  the  propeller  efficiency  increases  to  a  peak  and  then  decreases. 
Since  the  exact  propeller  rotational  speed  tested  by  01  et  al.  is  unclear,  it  can  only  be  compared 
to  the  trends  in  the  data.  The  present  data  agrees  well  with  that  shown  by  01  et  al.  At  a  rotational 
speed  of  n  =  8000  rpm,  the  propeller  was  tested  by  sweeping  the  advance  ratio  from  low  to  high 
values,  and  then  sweeping  from  high  to  low  values  to  examine  the  potential  for  hysteresis  in  the 
experiment.  As  can  be  seen,  there  is  not  a  noticeable  difference  between  these  two  sets  of  data. 
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Figure  49  Comparison  of  Present  Dynamic  Results  to  01  et  al.  (2008)  (APC  6.0  x  4.0) 
(a)  Coefficient  of  Thrust,  (b)  Coefficient  of  Torque,  (c)  Propeller  Efficiency. 
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5.6.4  Dynamic  test  results 

With  the  dynamic  results  validated,  data  was  collected  for  all  of  the  propellers.  Comparisons 
were  drawn  between  propellers  with  constant  diameter  and  varying  pitch  in  Figure  50  and 
between  propellers  with  constant  pitch  and  varying  diameter  in  Figure  51,  both  at  a  nominal 
rotational  speed  of  n  =  16000  rpm.  In  Figure  50,  propellers  with  larger  pitch  generally  had  larger 
coefficients  of  thrust  and  power,  and  the  windmill  state  occurred  at  higher  values  of  the  advance 
ratio,  which  indicates  that  larger  pitch  values  tend  to  allow  for  higher  airspeed.  Figure  50(c) 
shows  that  the  propeller  efficiency  decreases  with  increasing  pitch  for  lower  values  of  advance 
ratio,  and  the  peak  efficiency  occurs  at  higher  values  of  advance  ratio.  An  increase  in  pitch 
essentially  means  that  the  angle  of  attack  of  the  airfoil  is  higher,  which  should  increase  both 
thrust  and  torque  prior  to  reaching  stall.  In  Figure  51,  increasing  the  propeller  diameter  for  a 
given  pitch  tends  to  decrease  the  coefficient  of  thrust  and  the  coefficient  of  power,  and  the 
propeller  efficiency  increases  with  diameter  for  lower  values  of  advance  ratio.  Increasing  the 
diameter  for  a  given  rotational  speed  and  airspeed  actually  increases  the  thrust  and  torque  due  to 
the  increased  wingspan  of  the  propeller,  but  this  effect  is  negated  due  to  the  factor  of  D4  in  the 
denominator  of  Cj  and  the  factor  of  D5  in  the  denominator  of  the  Cp. 
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Figure  50  The  Effect  of  Varying  Propeller  Pitch  while  holding  Diameter  Constant  (Dynamic  Testing) 
(a)  Coefficient  of  Thrust,  (b)  Coefficient  of  Propeller  Power,  (c)  Propeller  Efficiency. 
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Figure  51  The  Effect  of  Varying  Propeller  Diameter  while  holding  Pitch  Constant  (Dynamic  Testing) 
(a)  Coefficient  of  Thrust,  (b)  Coefficient  of  Propeller  Power,  (c)  Propeller  Efficiency. 
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Figure  52  presents  results  for  square  propellers,  where  D  =  P,  from  01  et  al.  [95],  Selig  [96],  and 
the  present  experiment  for  a  fairly  wide  range  in  propeller  diameter  (4.0  <  D  <  18  inches).  01  et 
al.  had  conjectured  that  the  coefficient  of  thrust  should  collapse  for  square  propellers.  The  results 
are  grouped  from  small  to  large  propeller  diameter,  where  the  three  researchers  essentially 
covered  different  diameter  ranges.  The  data  does  not  appear  to  collapse  in  that  the  propeller 
diameter  still  affects  the  coefficient  of  thrust.  It  should  be  noted,  however,  that  this  set  of  data 
covers  multiple  manufacturers  and  propeller  types  (i.e.,  Thin  Electric,  Sport).  Geometric 
differences  in  the  blades,  such  as  the  twist  and  chord  distributions,  as  well  as  blade  stiffness  (and 
hence  deflection  under  loading)  may  significantly  affect  performance. 


Advance  Ratio 


Figure  52  Coefficient  of  Thrust  versus  Advance  Ratio  for  Square  Propellers  (D/P  =  1.0)  with  Propeller 

Diameter  Ranging  from  4.0  <  D  <  18  inches. 

Figure  52  presents  the  coefficient  of  thrust  for  the  same  family  of  propellers  (APC  Speed  400 
Electric).  It  was  assumed  that  using  a  group  of  similar  propellers  may  reduce  the  effects  of 
geometric  blade  variations.  In  this  case,  the  diameter-to-pitch  ratio  for  this  family  of  propellers 
has  a  relatively  small  range  (0.86  <D/P  <  1.5).  Figure  53(a)  shows  the  coefficient  of  thrust 
versus  advance  ratio  for  all  of  the  collected  data  for  the  APC  Speed  400  Electric  propellers  with 
an  uncertainty  level  of  ACt  <  20%.  As  can  be  seen,  the  results  are  not  correlated  well,  as 
witnessed  by  the  low  goodness  of  fit  parameter,  R2  =  0.529.  Figure  53(b)  shows  the  coefficient 
of  thrust  modified  by  DIP  plotted  against  advance  ratio,  where  the  goodness  of  fit  parameter  is 
improved  to  R2  =  0.680.  This  correlation  can  be  further  improved  by  modifying  the  advance  ratio 
by  DIP  as  seen  in  Figure  53(c),  where  the  goodness  of  fit  parameter  increases  to  R2  =  0.720, 
which  is  a  significant  improvement  over  the  original  data  correlation.  This  analysis  shows  that 
the  inclusion  of  the  propeller  pitch  into  the  correlations  for  the  performance  coefficients  can  be 
beneficial  for  this  family  of  propellers. 
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Figure  53  Coefficient  of  Thrust  vs.  Advance  Ratio  for  the  APC  Sport  400  Electric  Propellers  (ACt  <  20%) 
(a)  Original  Representation  of  Cf,  (b)  CT  Modified  by  D/P,  (c)  CT  and  J  Modified  by  D/P. 

Figure  54(a)  shows  the  coefficient  of  power  versus  advance  ratio  for  all  of  the  collected  data  with 
an  uncertainty  level  of  ACp  <  20%.  As  can  be  seen,  the  results  are  not  correlated  well,  as 


98 

Approved  for  public  release;  distribution  unlimited 


witnessed  by  the  very  low  goodness  of  fit  parameter,  R2  =  0.059.  Figure  54(b)  shows  a 
coefficient  of  power  that  is  modified  by  DIP,  which  improves  the  goodness  of  fit  parameter  to  R2 
=  0.538.  This  is  still  relatively  low,  and  probably  should  not  be  used  in  most  engineering 
analyses.  Figure  55  shows  a  similar  comparison  for  the  propeller  efficiency  versus  the  advance 
ratio.  Here,  the  diameter  to  pitch  ratio  was  used  to  modify  the  original  advance  ratio  to  increase 
the  goodness  of  fit  parameter  from  R2  =  0.938  to  0.983,  which  is  deemed  to  be  very  accurate  for 
most  applications,  especially  for  advance  ratios  less  than  J<  0.6. 
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Figure  54  Coefficient  of  Propeller  Power  versus  Advance  Ratio  for  the  APC  Sport  400  Electric  Propellers 

(ACP  <  20%) 

(a)  Original  Representation  of  CP;  (b)  CP  Modified  by  D/P. 
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Figure  55  Propeller  Efficiency  vs.  Advance  Ratio  for  the  APC  Sport  400  Electric  Propellers  (Ai/P<  20%) 
(a)  Original  Representation  of  r/P;  (b)  Advance  Ratio  Modified  by  D/P. 

5.7  Conclusions 

Twenty-four  propellers  in  the  range  of  4.0  <  D  <  6.0  inches  in  diameter  and  2.0  <  P  <5.5  inches 
in  pitch  were  tested  statically  and  dynamically  in  the  Wright  State  University  wind  tunnel  over  a 
wide  range  of  propeller  rotational  speeds  and  air  speeds.  A  detailed  experimental  procedure  for 
both  cases  was  employed  and  an  extensive  uncertainty  analysis  was  performed  on  the  resulting 
data.  The  experiments  were  validated  by  comparing  the  results  to  previous  works.  The 
repeatability  of  the  experimental  results  and  the  repeatability  of  the  manufacture  of  the  propellers 
were  proven  by  testing  three  duplicate  propellers  three  times  each.  Static  tests  were  performed  by 
varying  propeller  speed  from  n  =  4000  rpm  to  the  maximum  speed  limited  by  the  manufacturer’s 


101 

Approved  for  public  release;  distribution  unlimited 


specifications  or  the  maximum  motor  temperature.  Dynamic  tests  were  performed  by  holding  the 
propeller  speed  constant  and  varying  the  wind  tunnel  airspeed  and  thus  varying  the  advance  ratio. 

For  a  given  airspeed  and  rotational  speed,  the  thrust  and  torque  both  increased  with  propeller 
pitch  and  diameter,  as  expected.  For  static  testing,  the  coefficient  of  thrust,  coefficient  of 
propeller  power  and  the  total  propulsive  efficiency  increased  with  propeller  pitch  for  a  given 
propeller  diameter  at  all  rotational  speeds.  Alternately,  when  the  diameter  was  increased  and  the 
pitch  was  held  constant,  the  coefficient  of  thrust  and  coefficient  of  propeller  power  increased  and 
the  total  propulsive  efficiency  decreased.  Similar  results  were  found  for  dynamic  testing.  The 
square  (D/P)  propellers  that  were  tested  and  compared  to  the  results  in  the  open  literature  were 
not  correlated  well  with  propeller  diameter  alone,  which  was  possibly  due  to  variations  in  the 
twist  and  chord  distributions,  as  well  as  the  blade  stiffness.  Correlations  for  the  coefficient  of 
thrust,  propeller  power  and  propeller  efficiency  for  a  family  of  propellers  (APC  Speed  400 
Electric)  under  dynamic  testing  were  found  to  be  improved  when  the  original  coefficient  (and/or 
advance  ratio)  was  modified  by  D/P.  This  pennitted  the  pitch  to  be  accounted  for  in  the 
perfonnance  characteristics  of  propellers  in  the  range  of  0.86  <  D/P  <  1.5. 
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